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The formation of stars is a key process in astrophysics. Detailed knowledge of the physical mechanisms that 
govern stellar birth is a prerequisite for understanding the formation and evolution of our galactic home, the 
Milky Way. A theory of star formation is an essential part of any model for the origin of our solar system 
and of planets around other stars. Despite this pivotal importance, and despite many decades of research, our 
understanding of the processes that initiate and regulate star formation is still limited. 

Stars are bom in cold interstellar clouds of molecular hydrogen gas. Star formation in these clouds is governed 
by the complex interplay between the gravitational attraction in the gas and agents such as turbulence, magnetic 
fields, radiation and thermal pressure that resist compression. The competition between these processes deter- 
mines both the locations at which young stars form and how much mass they ultimately accrete. It plays out over 
many orders of magnitude in space and time, ranging from galactic to stellar scales. In addition, star formation is 
a highly stochastic process in which rare and hard-to-predict events, such as the formation of very massive stars 
and the resulting feedback, can play a dominant role in determining the evolution of a star-forming cloud. 
As a consequence of the wide range of scales and processes that control star formation, analytic models are 
usually restricted to highly idealized cases. These can yield insight, but the complexity of the problem means that 
they must be used in concert with large-scale numerical simulations. Here we summarize the state of modem star 
formation theory and review the recent advances in numerical simulation techniques. 



I. INTRODUCTION 

Stars are central to much of modem astronomy and astro- 
physics. They are the visible building blocks of the cosmic 
structures around us, and thus are essential for our under- 
standing of the universe and the physical processes that govern 
its evolution. At optical wavelengths almost all natural light 
we observe in the sky originates from stars. The Moon and 
the planets in our solar system reflect the light from our Sun, 
while virtually every other source of visible light further away 
is a star or collection of stars. Throughout the millenia, these 
objects have been the observational targets of traditional as- 
tronomy, and define the celestial landscape, the constellations. 
The most massive stars are very bright, they allow us to reach 
out to the far ends of the universe. For example, the most 
distant galaxies in the Hubble Ultra Deep Field are all charac- 
terized by vigorous high-mass star formation. Understanding 
the origin of stars, at present and at early times, therefore is a 
prerequisite to understanding cosmic history. 

Stars are also the primary source of chemical elements 
heavier than the hydrogen, helium, and lithium that were pro- 
duced in the Big Bang. The Earth, for example, consists 
mostly of iron (32%), oxygen (30%), silicon (15%), magne- 
sium (14%) and other heavy elements (Mo rgan and Anders) 
1980| l. These are produced by nuclear fusion in the interior of 
stars, and enriching gas to the chemical composition observed 
today in our solar system must have required many cycles of 
stellar birth and death. 



Today we also know that many stars harbor planetary sys- 



tems around them, about 300 are known as of fall 2008 ( Udry 



and Santos 2007). The build-up of planets is intimately cou- 
pled to the formation of their host stars. Understanding the 
origin of our solar system and of planets around other stars 
has a profound impact on how we see our position in the uni- 
verse. Questions whether we are alone, or whether there is life 
elsewhere in the cosmos are of broad interest to all of us. 

Stars and the planetary systems they may harbor are born 
in turbulent interstellar clouds of molecular hydrogen with 
a small fraction of dust mixed in. At optical wavelengths, 
we see these clouds as dark patches of obscuration along the 
band of the Milky Way. The dust component blocks the light 
from stars further away. At far-infrared, sub-millimeter, and 
radio wavelengths, however, the dust becomes increasingly 
transparent and we can look into these clouds. These ob- 
servations reveal extremely complex morphological and kine- 
matic structure, where patches of cold high-density gas are 
interspersed between regions of low-density warmer mate- 
rial ( |Ferriere| |2001| l. It is thought that this complicated tex- 
ture is caused by supersonic turbulence that is generated by 
large-scale gravitational motions in the galaxy (such as spi- 
ral density waves) or by energy and momentum input from 
stars themselves (Elmegreen and Scalo ' 2004{ [Mac Low and] 
pGessen, ,2004, .Scalo and Elmegreen 2004[ ). 

Within molecular clouds, supersonic turbulence and ther- 
mal instability lead to a transient, clumpy structure. Some of 
the resulting density fluctuations exceed the critical mass and 
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FIG. 1 Star forming region NGC 602 in the Large Magellanic Cloud, 
observed at optical and infrared wavelenths. The intense radiation 
from high-mass stars in the center of the young cluster has carved 
a cavity into the surrounding parental molecular cloud. Elephant 
trunk-like dust pillars that point towards the hot blue stars are the 
signs of this eroding effect. (Image courtesy of NASA, ESA, and the 
Hubble Heritage Team ) 



density of gravitational stability. These clumps begin to col- 
lapse, their central density increases rapidly, and eventually 
they give birth to new protostars. Clusters of stars form in 
large regions that become unstable, within which contraction 
involves multiple collapsing cloud cores. A number of re- 
cent reviews have discussed various aspects of stellar birth in 



these clouds ( Ballesteros-Paredes et al.\ 


20071 Larson 2005| 


Mac Low and Klessen||2004| McKee and 


Ostriker 2007 Zin-| 


necker and Yorke 2007|l. 



Stars and their parental clouds are connected via a number 
of feedback loops. Stars of all ages radiate and will thus heat 
up the gas in their vicinity. By doing so they influence sub- 
sequent star formation. Massive stars emit photons at ultravi- 
olet wavelengths, creating bubbles of hot, ionized gas around 
them (Be tidier ef"ar||2007[|Hoare et ar||2007| l, as illustrated 
in Figure [T] These so-called Hll regions are likely to quench 
further star formation in their interior, and thus set the star- 
formation efficiency in the region. The collective action of 
many Hll regions can destroy entire molecular clouds, and 
thus have the potential to influence the star-forming proper- 
ties of galaxies on larger scales. The combined effect of large 
numbers of supernova explosions is another important mech- 
anism for driving the supersonic turbulence ubiquitously ob- 
served in the galactic gas. By the same token, however, these 
feedback processes may trigger the birth of new stars. The 
very same processes that terminate star formation in one loca- 
tion may compress gas somewhere else in the galaxy, leading 



to new star formation. 

The density contrast between typical cloud densities and 
the hydrogen-burning centers of the final stars is enormous, 
about 24 orders of magnitude, and so is the coiTesponding 
spatial range (roughly 8 orders of magnitude). In addition to 
the large dynamical range, many different physical processes 
play a role at the various stages of the contraction process. On 
global scales we need to describe the formation of molecular 
clouds via large-scale flows of mostly atomic gas in a galac- 
tic disk. Internal turbulent compression in the star-forming 
cloud then sets the initial stage for the protostellar collapse 
of individual objects. The thermodynamics of the gas, and 
thus its ability to respond to external compression and conse- 
quently to go into collapse, depends on the balance between 
heating and cooling processes. Magnetic fields and radiative 
processes also play an important role. Modeling star forma- 
tion adequately therefore requires the accurate and simultane- 
ous treatment of many different physical processes over many 
different scales. 

In the past, progress has only been achievable by divid- 
ing the problem into smaller bits and pieces and by focus- 
ing on few physical processes or single scales only. Today, 
however, algorithmic advances and increasing computational 
power permit a more integrated approach to star formation. 
For the first time we are able to combine, for example, magne- 
tohydrodynamics with chemical and radiative processes, and 
apply these numerical schemes to real astrophysical problems. 
It is this integrated view of stellar birth that is at the heart 
of this review. Our goal is to provide an overview of the re- 
cent advances in star formation theory with a special focus on 
the numerical aspects of the problem. We do not aim to be 
complete, for this we refer the reader to the recent reviews in 
this fi eld (|Mac Low and Kless enl|2004[ |McKee and Ostriker[ 
|2007[ [Zinnecker and Yorke| |2007| l. Instead we will focus on 
three selected topics where we think numerical studies have 
had the largest impact and where we think our understanding 
of the physical processes that initiate and regulate stellar birth 
evolves most rapidly. We begin with the large scales and dis- 
cuss the formation of molecular clouds in galactic disks and 
the numerical requirements and methodologies needed to do 
so consistently in Section [ll] We then zoom into individual 
star-forming regions and examine the transition from cloud 



cores to stars in Section III As a third focus point, we discuss 
in Section|lV]the effects of stellar feedback and examine how 
it alters the star formation process. Finally, in Section [V] we 
summarize and speculate about the future of numerical star 
formation research. 



II. THE SITES OF STAR FORMATION: MOLECULAR 
CLOUDS 

A. Phenomenology of Molecular Clouds 

In regions of the interstellar medium (ISM) that are suf- 
ficiently dense and well-shielded against the dissociating ef- 
fects of interstellar ultraviolet radiation, hydrogen atoms bind 
to form molecules. Star formation appears to occur exclu- 
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FIG. 2 Molecular cloud complex in the constellation Perseus. The 
image shows the distribution of CO line emission at radio wave- 
lengths. This is a good tracer of total gas mass. Clearly visible is 
the highly complex and filamentary morphological structure of the 
cloud. (Image from \Sun et aL] ) |2006p .) 



sively within this molecular phase of the ISM. Molecular hy- 
drogen is a homonuclear molecule, so its dipole moment van- 
ishes and it radiates extremely weakly. Direct detection of 
cold interstellar H2 is generally possibly only through ultravi- 
olet absorption studies, such as those made by the the Coper- 



nicus (Spitzer and Jenkins 19751 and Far Ultraviolet Spec- 
troscopic Explorer (FUSE) (Moos ef^r||2000 ; Sahno w et al.\ 
|2000| l satellites. Due to atmospheric opacity these studies can 
only be done from space, and are limited to pencil-beam mea- 
surements of the absorption of light from bright stars or active 
galactic nuclei. Note that rotational and rovibrational emis- 
sion lines from H2 have also been detected in the infrared, 
both in the Milky Way and in other galaxies. However this 
emission comes from gas that has been strongly heated by 
shocks or radiation, and it traces only a small fraction of the 
total H2 mass (e.g. |van der Werf[|2000[ ). Due to these limita- 
tions, the most common tool for study of the molecular ISM 
is radio and sub-millimeter emission either from dust grains 
or from other molecules that tend to be found in the same lo- 
cations as H2. The most prominent of these is CO, although 
other tracers such as HCN are beginning to come into wide 
use. 

As of this writing, the Milky Way and a few dozen Lo- 
cal Group galaxies have been mapped in the J = 1 ^ 
or 2 ^ 1 rotational transitions of CO at a resolution bet- 
ter than 1 kpc ( |Bigiel et aI]j 20Q8 , Blit z et a/.|[20 07 ; Bolatto 

i87hEng argiola et a/. ','2003 ; Fukui 
T999; Leroy ef a/(|2008; Rosolowsky and BHtz, 



et al. 



et al. 



2005 



[2007| p 
p008][ 



Rosolow sky et a/.[ |2003|]Solomon et al] |1987| |WaFl 



ter et al., ,2008| ), and a large number of more distant galax- 
ies have been imaged at lower spatial resolution. The frac- 
tions of the ISM within the molecular phase in these galax- 
ies ranges from no more than a few percent in low-surface 
density dwarfs to near unity in giant high-surface density sys- 
tems. The highest molecular fractions are generally found in 
the parts of galactic disks with the highest total gas surface 



densities, but in the most actively star-forming galaxies the 
molecular fraction can reach ~ 90% even integrated over the 
entire galaxy ( |Iono et al.\ |2005| [Mirabel and Sandersj [T989| . 
In all of the nearby galaxies where high resolution observa- 
tions are possible the molecular gas is largely organized into 
giant clouds (the so called giant molecular clouds or GMCs) 
of mass w 10* — 10^ Mq with average densities ^ 100 H2 
molecules per cm'^, separated by a more diffuse intercloud 
medium. In the Milky Way and galaxies of lower density this 
medium is mostly atomic or ionized hydrogen, while in the 
densest nearby galaxies, such as M64 ( [Rosolowsky and BIitz[ 
[2005[ ), it is also molecular. 

Molecular clouds across the Local Group all seem to dis- 
play a number of properties in common. First, when stud- 
ied with high spatial resolution clouds, they exhibit extremely 
complex and often filamentary structure, with column densi- 
ties and corresponding 3-D densities that vary by many or- 
ders of magnitude (see Figure |2] or Table |l]l. Nevertheless, 
when observed with low resolution, to within factors of a 
few all molecular clouds seem to have a similar mean sur- 
face d ensity of ~ 100 Mp, pc~^ corresponding to 0.035 g 



cm"^ (Bolatto et al. 2008 [Heyer et al. 2008 1. The con 



stant surface density of molecular clouds is known as one 
of Larson's Laws (Larson 1981[ l, although there are a num- 
ber of caveats with these relations and their interpretation 
( [Ballesteros-Paredes and Mac Low|[2002[ ). Second, the clouds 
all display linewidths much greater than would be expected 
from thermal motion, given their inferred temperatures of 
10 — 20 K. The observed linewidth is related to the size of 
the cloud by 



(TiD = 0.5 



L 



1.0 pc 



0.5 



km s 



(1) 



where ctid is the one-dimensional cloud velocity dispersion 



and L is the cloud size (Bolatto ef a/. 2008 Heyer and Brunt 



[2004[ [Solomon et a/.||1987 ). This is another one of Larson's 
Laws. These non-thermal linewidths have been interpreted as 
indicating the presence of supersonic turbulent motion, since 
both the low observed star formation rate (see below) and the 
absence of inverse P-Cygni line profiles indicates that they are 
not due to large-scale collapse. If one adopts this interpreta- 
tion, then from these two observed relations one can directly 
deduce the third of Larson's Laws, which is that giant molecu- 
lar clouds have virial parameters (Bertoldi and McKee 1992| l 



GM 



1, 



(2) 



where M is the cloud mass. This indicates that these clouds 
are marginally gravitationally bound, but with enough inter- 
nal turbulence to at least temporarily prevent global collapse; 
whether they are truly in virial equilibrium is a topic that we 
discuss in detail below. (There is also a population of molec- 
ular clouds with virial ratios ^ 1, but these have masses 
<C lO"' Mq, and do not appear to host star formation (Heyer| 
[ef a/.|[200T] ).) These relations appear to partially or fully break 
down in starburst galaxies with very high surface densities, 
where for example the molecular gas velocity dispersion can 
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TABLE I Physical properties of molecular cloud and cores" 





molecular 


cluster- 


protostellar 




cloud 


forming 


cores 






clumps 




Size (pc) 


2-20 


0.1 - 2 


<0.1 


Density (n(H2)/cm^) 


10^ - 10* 


10^ - 10^ 


> 10" 


Mass (M0) 


10^ - 10* 


10- 10^ 


0.1 - 10 


Temperature (K) 


10-30 


10-20 


7-12 


Line width (kms^^) 


1 - 10 


0.3-3 


0.2-0.5 


Column density 








(g cm"^) 


0.03 


0.03- 1.0 


0.3-3 


Crossing time (Myr) 


2-10 


<1 


0.1 - 0.5 


Free-fall time (Myr) 


0.3-3 


0.1 - 1 


<0.1 


Examples 


Taurus, 


L1641, 


B68, L1544 




Ophiuchus 


LI 709 





reach ~ 100 km s^^ ( jPownes and Solomon 1998 i, but it is 
unknown whether there are analogous relations under these 
higher density conditions. 



The presence of supersonic turbulence in approximate virial 
balance with self-gravity indicates that in molecular clouds 
the turbulent and gravitational energy densities are of the same 
order of magnitude, and both greatly exceed the thermal en- 
ergy density. If molecular clouds form in large-scale con- 



vergent flows (as we argue below in Section II. B 1, then sur- 
face terms from ram pressure can also be significant and need 
to be considered in the virial equations ( jBallesteros-Paredes} 
2006| l. There is also one further energy reservoir that we must 
mention, magnetic fields. The gas in molecular clouds is a 
weakly ionized plasma that is tied to magnetic field lines. Ob- 
servations using Zeeman splitting ( [Crutcher] |1999[ [Troland 
and Crutcher, 2008) and the Chandrasekhar-Fermi effect ( Lai 
eTa/.l |2001| |2002| ) indicate that the field strength lies in the 
range from a few to a few tens of fj,G. The exact value varies 
from region to region, but in general the magnetic energy den- 
sity appears comparable to the gravitational and turbulent en- 
ergy densities. One can describe this state of affairs in terms 
of magnetic criticality. If the magnetic field threading a cloud 
is sufficiently strong, then it cannot undergo gravitational col- 
lapse no matter what external pressure is applied to it, as long 
as it is governed by ideal magnetohydrodynamics (MHD). A 
cloud in this state is called subcritical. In contrast, a weaker 
magnetic field can delay collapse, but can never prevent it, 
and a cloud with such a weak field is called supercritical 



(Mouschovias and Spitzer 1976 1. Observations indicate the 
molecular clouds are close to being, but not quite, magneti- 
cally subcritical (' Crutcher et al.\ |2008a| l. For further discus- 
sions, see Section pI.A.2| 



B. Cloud Timescales and Cloud Formation 

1 . Characteristic Timescales for IVIolecular Clouds 

We can learn a great deal about molecular clouds by con- 
sidering the timescales that govern their behavior. Because 
molecular clouds span a huge range of size and density scales, 
and their evolution times reflect this range, it is convenient to 
normalize all discussion of timescales to the free-fall time, de- 
fined as the time that a pressureless sphere of gas with some 
initial starting mass density p requires to collapse to infinite 
density under its own gravity: tff — 3tt / (32G p) . For a 
cloud for which the virial parameter avir ~ 1, this is roughly 
half the cloud crossing time (Tan et al. 20061, defined as 



the ratio of the characteristic size to the velocity dispersion 
tcr = i/ciD- The timescales tg or tcr define the characteris- 
tic timescales on which behavior driven by gravity or limited 
by the internal gas signal speed can operate. For a molecular 
cloud detected via CO emission, with a mean number den- 
sity n « 100 cm"-^, tff w 3 X 10^ yr. We now define some 
other useful timescales that can be determined from observa- 
tions, and which yield strong constraints on how molecular 
clouds must behave. Any complete theory of star formation 
in molecular clouds must be able to explain each of the three 
timescales we describe. 

The Gas Depletion Time. Perhaps the most fundamental ob- 
servational timescale for molecular clouds is the gas depletion 
time <dcp, which is defined as the ratio of the mass of a molec- 
ular cloud (or population of clouds) to the star formation rate. 
This defines the time that would be required to convert the 
cloud completely into stars at its observed star formation rate 
(SFR), assuming this rate is constant over time. Estimating 
this number immediately yields a striking conclusion, which 
is perhaps the most basic problem in star formation. The disk 
of the Milky Way contains « 10^ Mq of molecular gas, and 
the observed star formation rate is only a few Mo/year, so 



the gas depletion time tdep must be a few hundred Myr ( Mc 



Kee 1999 Zuckerman and Evans 1974 1, roughly 100 times 
the free-fall time. (Note, that if we compare this timescale 
with the age of the Milky Way of close to 10^° yr, we con- 
clude that a continuous inflow of fresh gas is required if the 
cuiTent SFR is at all representative and if we assume we are 
not living in times when the Milky Way is running out of gas 
soon. This problem gets worse if we consider proposals that 
the SFR was higher in the past ( |Madau et al] |1998| ).) One 
can repeat this exercise for populations of molecular clouds in 
both the Milky Way and in other galaxies ( Bigiel ef a/. , 2008 ), 
using a variety of indicators of the star formation rate ( |Kenni-| 
|cutt| |1998| l, and using a similarly wide variety of techniques 
to estimate masses of molecular clouds with various densities. 
Doing so yields the data shown in Figure |3] In this Figure the 
X-axis indicates the characteristic density to which a particu- 
lar method of measuring molecular gas is sensitive, and the 
y-axis shows the ratio tff/idop for that gas. The fact that this 
ratio is w 1 % for low density gas and either remains constant 
or slowly increases to at most 10% at high densities indicates 
that the conversion of gas into stars must be inefficient or slow, 
in the sense that no more than a few percent of the total mass 
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in molecular clouds in a galaxy can be converted into stars 



per free-fall time ( Krumholz et al. 2006 Krumholz and Tan 
|2007| l. This discrepancy is at the heart of any star formation 
theory, but before we can address it we must consider some 
other important timescales. 
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FIG. 3 Ratio of molecular cloud free-fall times to depletion times, 
versus the characteristic hydrogen density n to which the indicated 
tracer of the molecular gas is sensitive. The depletion time is defined 
as the time that would be required to convert all of the gas into stars 
at the observed star formation rate. Each data point represents a sur- 
vey of molecular clouds using a different tracer that probes gas of 
different densities, from CO (1^0) emission to CS (5^4) emission, 
which yields only an upper limit. IRDCs stands for infrared dark 
clouds, which are detected in infrared absorption, and ONC stands 
for the Orion Nebula Cluster, a single nearby star-forming region. 
(Adapted from \K'rumholz and Tan\li2007\ .) 



The Molecular Cloud Lifetime. The gas depletion time tells 
us how long it would take to convert a molecular cloud into 
stars completely. The actual lifetime of the cloud, how- 
ever, is considerably shorter. Most of the cloud's mass is never 
converted into stars. Instead it participates in the perpetual cy- 
cle that connects the molecular, atomic, and ionized phases of 
the ISM ( Ferriere[ |2001| l. Molecular clouds form out of the 
atomic gas as discussed below (Section II. B. 2 1, convert some 
of their mass into stars, and then dissolve either by internal 
feedback (Section [TVl or large-scale dynamics. The total du- 
ration of this process is very hard to estimate. 

In external galaxies, estimates of molecular cloud lifetimes 
are usually obtained from statistical relations between the lo- 
cation of molecular clouds and young star clusters. Star clus- 
ters can be age-dated, and determining at what ages they cease 
to be located preferentially close to molecular clouds gives an 
estimate of how long molecular clouds live once they form 
visible star clusters. Correcting for the population of molecu- 
lar clouds that have not yet produced optically-visible clusters 
then gives an estimate of fufo. In the LMC this is tiifc ~ 27 
Myr ([Blitz et a/.|[2007| and in M33 it is tnfe ~ 20 Myr ( [Ei>] 
Igargiolaef a/.[ 2003 1, so <iifo ^7—10 ts, with a factor of ~ 2 
uncertainty. Due to the sensitivity limits of the observations, 
these estimates apply only to GMCs ^ 10^ Mq or larger. 

In star-forming regions within ~ 500 pc of the Sun, we can 
obtain ages estimates using individual young stars (as opposed 
to star clusters). Since stellar populations older than about 



5 X 10^ yr are generally not associated with molecular gas 
anymore, this technique suggests m olecular cloud lifetimes 
substantially shorter than 10*^ years ([Ballesteros-Paredes and 



Hartmann 2007 Hartmann et al. 2001 1. Placing stars on 



the Hertzsprung-Russell diagram results in stellar age spread 



from IMyr up to 3 Myr (Hartmann 2003; Huff and Stabler 



|2006t |Palla and Stabler, ,2000j , with considerable uncertainty 
but consistent with the above number. Unfortunately these 
nearby clouds have all masses well below ^10^ Mq, so there 
is essentially no overlap between this population and the ex- 
tragalactic one. Probably in part as a result of this selec- 
tion effect, these regions are also denser than the giant clouds 
we can observe in external galaxies, and consequently have 
lower free-fall times, so tnfe ~ 1 — 3 Myr corresponds to 
iiife ~ 1 - lOiff dTanef a/.|[2006l ). 

The Lag Time. The third timescale describing molecular 
clouds is the lag between when they form and when they be- 
gin to form stars, which we call the lag time iiag. For molecu- 
lar clouds in the solar neighborhood (out to 800 pc), the ratio 
of star-forming clouds over those without clear signs of star 
formation ranges between 7 and 14 ( [Ballesteros-Paredes and| 
~" ' 2007[). Together with a median age of the young 



Hartmann 



stars in these regions of 1 — 2 Myr (Hartmann '2003 ; Pallaj 
|and Stahler| 12002 ) this entails a lag between cloud formation 
and onset of star formation of at most tiag ~ 1 Myr, i.e. stars 
begin to form immediately after (or even during) the forma- 
tion of the parental cloud. This is consistent with some extra- 
galactic observational evidence that the spatial gap between 
spiral shocks in H I gas and bright 24 fjm emission down- 
stream of the shock, presumably tracing star formation, corre- 



sponds to a lag time tiag = 1—4 Myr (Tamburro et al. 2007 1. 



Before proceeding based on this conclusion, however, it is 



worth mentioning a caution. In both M33 (Engargiola et al. 



|2003[ l and the LMC ( Blitz ef al.\\2001\ , the ratio of molecular 
clouds that have associated H II regions (detected either via 
Ha or radio continuum emission) to those that do not is sig- 
nificantly smaller than the local ratio of star-forming clouds 
to non-star-forming ones: 3 — 4 instead of 7 — 14. This im- 
plies lag times of ^ 7 Myr, roughly 2 — 3iff, between GMC 
formation and H II region appearance. While the extragalac- 
tic clouds used for this measurement are much larger than the 
solar neighborhood ones and have free-fall times of ~ 3 Myr 
instead of ~ 1 Myr, they are comparable in size to the clouds 



probed by the geometric technique ( Tamburro et al. 2007 1. 
The discrepancy between ty^g — 2 — 3tff and tiag^ l^ff be- 
tween these techniques is therefore real. Its origin is unclear. 
One possibility that it is simply a result of the different criteria 
used to measure the onset of star formation: infrared emission 
versus H II regions. This possibility has yet to be quantita- 
tively evaluated, however In the absence of an explanation 
of the discrepancy, we tentatively conclude based on the IR 
data that tiag is indeed short. Because of this extremely short 
timelag, any property the molecular cloud has to allow star 
formation, i.e. the strong density variations including small 
filling factors, and the supersonic turbulence, must come from 
the formation process of the cloud itself 

Strong density variations within molecular clouds are not 
only observed (Section III. A. 1 Table |l]l, they also are phys- 
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ically mandated. The free-fall time for a spherical cloud of 
uniform density does not depend on the radius. Thus if one 
neglects pressure forces, material at the edge of the cloud 
would arrive at the same time at the center as material close 
to the center, making it virtually impossible to form isolated 
stars. "Distributed" star formation can only occur if the cloud 
acquires high, localized density seeds during its formation 
process or similarly if pre-existing density fluctuations exist 
that become strongly amplified, such that the local contrac- 
tion time is substantially smaller than the global one ( ,Burkert 
and Hartmaiinl |2004[ [Heitsch and Hartmann] |2008t |Heitsch 
et al. 2008b| l. The distribution of these high-density regions 
determines the degree of clustering of the resulting stars, with 
over-densities that are correlated on small length scales lead- 
ing to isolated stars or small clusters, while over-densities cor- 



related on large scales produce large clusters (Klessenl 2001a 



Klessen et a/. | [2000). One possible explanation of the pres- 
ence of small scale density inhomogeneity in newborn molec- 
ular clouds is that they form from atomic gas that is thermally 



bistable (Section II. B. 2 1. The onset of thermal instability leads 
to a wide-spread distribution of small-scale non-linear den- 
sity fluctuations on very short timescales. This could explain 
how strong density variations occur on small scales even if 
molecular cloud turbulence is driven on large scales by the as- 
sembly of the cloud in the Galactic disk. In addition, GMC's 
are highly filamentary and show sheet-like morphology ( [Blitz | 
et al. |2007| l. This means that also boundary effects are im- 
portant and any pre-existing initial fluctuations are more eas- 
ily amplified compared to models mentioned above that as- 
sume spherical symmetry ( [Burkert and Hartmann[|2004[|Hart-| 
mann and Burkertl|2007| l. We note, however, that the smaller 
structures in their interior that form star clusters do appear to 
be more centrally concentrated ( [Mueller et al.\ \2QQ2\ [Plume[ 
^Far][T997l [Shirley et a/.|[2003) . 



2. Molecular Cloud Formation 

The strict observational limits on iiag, the time between 
when molecules first appear and when star formation begins, 
have recently led to a focus on the process of molecular cloud 
formation. This can be split into three issues, namely the ac- 
cumulation problem, the chemistry problem, and the issue of 
rapid fragmentation. We will discuss each of them in turn. 

The Accumulation Problem. Building a molecular cloud 
requires assembling a column density high enough for the 
dust to shield the cloud against the ambient UV-radiation, 
and thus to allow CO formation. (H2 forms earlier, because 
it self-shields very efficiently ("Bergin et al.', '2004' 'Draine 
and Bert oldi) [T9961 [Krumholz et al. , 2008 ; van Dishoeck and 
Black[[1986|l. However, we "see" the cloud only once it con 



tains CO.) Dust-shielding becomes efficient at column densi- 
ties of approximately N = 2 x 10^^ cm. If we were to accu- 
mulate our prospective cloud at flow parameters typical for the 
inter-arm ISM, namely at densities 1 cm^-^ and velocities of 
lOkms^^, it would take ^60Myr to accumulate the shield- 
ing column density. This timescale is too large, given the max- 



However, this seemingly compelling argument is not ap- 
plicable, since it only addresses a one-dimensional situation, 
resulting in a sky-filling molecular cloud. Various three- 
dimensional solutions have been suggested to allow molecular 
cloud formation within realistically short times. Probably the 
oldest relies on the Parker instability ( [Parker[ [l966| l, a mech- 
anism describing the buyoant rise of evactuated parts of flux 
tubes in a stratified Galactic disk. The rising flux tube leads 
to material falling into the valleys, thus accumulating molec- 
ular clouds. In combination with a magnetohydrodynamical 
Rayleigh-Taylor instability to trigger the initial evacuation in 
spiral shocks fMouschovias et al. 19741), this scenario pre- 
dicts accumulation timescales of 30 Myr and typical cloud dis- 
tances of 1 kpc along spiral arms. However, more recent nu- 
merical simulations of magnetized galactic disk gas dynamics 
( Dobbs and Price , 2008 ; Kim and Ostriker 2006| l only found 
weak signatures of a Parker instability acting along a spi- 
ral arm. Instead they identified the Magneto-Jeans instability 



(Eknegreen 1987 Kim and Ostriker 2001 1 as a more domi- 
nant accumulation mode in a magnetized disk. This instability 
is driven by in-plane magnetic fields countering the stabilizing 
effects of the Coriolis force, allowing self-gravitating contrac- 
tions of overdense regions. Although we know that magnetic 
fields exist in galactic disks, it is also possible for molecular 
clouds to form on reasonable timescales in non-magnetized 
models provided the disk is globally gravitationally unstable 
([Dobbs et al. . 2 008 ; Li et al. 2005 2006 ; Tas keref a/.| [20081 
[Tasker and Tan|[2009p . 

On a more local scale, the interstellar medium is filled with 
flows of various types, many of which result in piling up ma- 
terial through shocks. While it is not always easy to iden- 
tify specific driving sources in all cases (see, however. 



'Nigral 



et al. 2008 ; Pate l et aL\ [1998 )), this is not surprising given 



the complexity expected as a result of the interaction of neigh- 
boring flows. In addition, the presence of massive molecular 
clouds well out of the Galactic plane (e.g., Orion) clearly sug- 
gests the need for some kind of driving. Given the extensive 
impact that massive stars have on the interstellar medium - 
Hll regions, stellar wind impacts, and ultimately supernova 
explosions - it is difficult to see how further creation of new 
star-forming locales by flows with scales of several to tens of 
pc (or even kpc, in the case of spiral arms) could be avoided. 
The notion of expanding shells piling up material has led to 
a "collect & collapse" model (Elmegreen 1998 ), connecting 
the formation of molecular clouds to energetic events in the 
ISM. 

The Chemical Timescale Problem. The formation of molec- 
ular hydrogen on dust grains - the main branch under Galactic 
conditions - is limited by three factors, namely the shield- 
ing from dissociating UV radiation, the dust temperature and 
the gas density ( [Tielens[ [2005 ). In the extreme case of zero 
H2 abundance in the assembling flows, dust shielding column 
densities need to be built up. However, due to its abundance of 
lines, H2 strongly self-shie lds already at column densities of 
N « 10^^ cm-^ ([Draine and Bertoldi ,1996, ; Tielens , ,2005j 



imum lifetimes of GMCs of < 30 Myr (Section II.B.l 



[van Dishoeck and Black 1986| l. Thus, alread small traces of 
H2 in the accumulating flows will help to lower the timescales 
for molecule formation (Pringle et al. 200I| l. The dust tem- 
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perature determines the efficiency of H2 formation on dust 
grains (i.e. what fraction of the accreted hydrogen atoms react 
to form H2). The efficiency is of order unity for Tdust < 25 K, 
but drops sharply above that, although it may remain as high 
as ~ 0.1 even up to Tdust = 1000 K (Cazaux and Tielens 



[2004 ). The timescale to reach equilibrium between formation 
and dissociation is given by w 10^ /n yr (HoUenbach et al. 
[T97T] l. 

Because of the complexity of the chemical reaction net- 
works, most cloud chemistry models are restricted to one- 
dimensional geometries, albeit in different environments such 



as molecule formation behind shock fronts (Bergin et al. 



2004^ or in (close to) quiescent clouds (Goldsmith et al. 



2004J 
2007^ 



20071 this is a good approximation. Models including hydro- 



dynamical effects such as shock compressions (Bergi n et al.\ 
[2004 ) or turbulence ( Glover and Mac Low 2007b ) predict 



shorter formation timescales - on the order of a few 10^ years 
- than static models ( [Glover and Mac Low[ [2007a[ ), empha- 
sizing the role of density variations and turbulent flows for the 
chemistry of the interstellar medium. 

The Fragmentation Problem. The key to the observation- 
ally mandated rapid onset of star formation is to provide the 
parental cloud with high-amplitude (non-linear) density per- 
turbations during its formation. We will describe a scenario 
of flow-driven cloud formation and rapid star formation in the 
context of cloud formation in spiral arms. While differing in 
details, other environments such as cloud formation in galaxy 
mergers or in the Galactic molecular ring, are subject to simi- 
lar physical constraints and may work along similar lines, al- 
though this is an issue still to be explored. 

Rapid fragmentation of the accumulating flows results from 
the combined action of heating and cooling processes in the 
ISM. In the diffuse, warm interstellar gas, at densities of 
n K, 1 cvcT^ photo-electric heating by dust grains dominates, 
while in denser, higher extinction gas, heating by cosmic rays 
is more significant, becoming dominant deep in the interior of 
molecular clouds. In the regime dominated by photo-electric 
heating, the total heating rate per hydrogen nucleus varies by 
at most a fact or of 10 over a range of de nsities from n w 1 
to 10"^ cm~3 ( |Wolfire et al.\ [1995[ [2003| . Cooling rates be- 
low 10"* K depend strongly on the abundances of heavy ele- 
ments, but have only a weak dependence on temperature at 



T > 100 K (Dalgarno and McCray 19721. Moreover, in the 



warm, diffuse ISM, the energy radiated in the dominant cool- 
ing lines, such as the hyperfine-structure line of singly ion- 
ized carbon at 158 /im, scales as the -n?, while the heating rate 
depends (roughly) on n. Starting from an equilibrium situa- 
tion, a small density increase thus leads to a cooling instability 
( [Field[ [T965| l. If the size of the density perturbation is small, 
with a sound-crossing time that is less than its cooling time, 
then as the gas cools, its density will increase owing to com- 
pression from the surrounding warmer medium. If the temper- 
ature dependence of the cooling rate is weak, then the increase 
in the cooling rate produced by the growing density is greater 
than the decrease caused by the falling temperature. And so 
the perturbation cools with ever faster growing density. This 
process will stop when the density dependence of the cooling 
rate changes, for instance, if the level populations of the dom- 



inant coolants reach their local thermodynamic equilibrium 
values. In this case the cooling rate scales only as n. It will 
also stop when the temperature dependence of the cooling rate 
becomes steeper, as will naturally occur at low temperatures. 
In the ISM, both effects are important, and the therm al insta- 
bility vanishes once n ^ 100 cm^'' and T 50 K ( [Wolfire 
\et al.\ [2003[ ), resulting in a two-phase structure of interstellar 
gas, with a warm diffuse phase occupying large volumes, and 
a cold, dense phase with small filling factors in rough pres- 
sure equilibrium ( [Burkert and Lin[ [2000; ,Heiles and Troland[ 
[2003] ). 

This thermal instability (in its various forms) is at the heart 
of the rapid fragmentation of the accumulating flows. It has 
been studied in various contexts, such as in generally turbulent 
media ([Audit and Hennebelie '2005', 'Heitsch et al.\ '2006b 
Kritsuk and Norman, 2002a b; Robertson and Kravtsov 



2008[), in cloud formation behind shockfronts ( |Koyama anc 



Inutsuka[ [2002[ [2000) , or in collisions of gas streams in spi 



ral arm shocks or driven by e.g. expanding supernova shells 
(Audit and Hennebelle, 2005; Heitsch et [2008bj [Hei>l 
inebelle et al., ,2008, ,Vazquez-Semadeni et al.\ 20071). The 
strength of the thermal instability derives from a combination 
with dynamical instabilities breaking up coherent shock fronts 



and shear-flow instabilit ies (Heitsch et al. 2006b Vazquez 
[Semadeni eFaL] [20061 [Vishniac[ [1994^ and from the fact 
that its growth timescales are substantially shorter than those 
of the hydromagnetic and gravitational instabilities involved 
( [Heitsch efaT|[2008a| l. 

The principal effect of this rapid thermal fragmentation is 
best gleaned from the evolution of the free-fall times in the 
forming cloud ( [Heitsch and Hartmann[|20(j8] l. Figure [4[ shows 
the distribution of free-fall times against evolution time in a 
molecular cloud being formed by two colliding flows of warm 
neutral hydrogen gas. Initially, the bulk of the cloud mass is at 
free-fall times longer than the simulation duration. At around 
7 Myr, a substantial mass fraction of the cloud has dropped to 
free-fall times as short as 3 Myr Substantial CO has formed 
by 10 Myr, while star formation sets in at « 11 Myr, when 
noticeable mass fractions appear at free-fall times substan- 
tially shorter than those in the bulk of the cloud. 

One of the key realizations in contrast to earlier models of 
turbulent fragmentation using periodic boxes is that the finite 
cloud geometry is crucial not only for the rapid onset of star 



formation (Heitsch et al.\ [2008b[ [Vazquez-Semadeni et al. 



2007 1, but also for the rapid formation of CO to overcome 



the otherwise stringent limitations set by the inflow density 
and speed ( Heitsch and Hartmann 2008 ). Thus these models 



bridge a gap between the large-scale simulations of Galactic 
disk dynamics discussed earlier, and the detailed models of 
turbulent fragmentation to be discussed in Sec [HI] The large- 
scale models do not have sufficient resolution to address the 
fragmentation and internal dynamics of the resulting molecu- 
lar clouds, while models on smaller scales generally have to 
make simplifying assumptions about the boundary conditions. 

What's Missing in Cloud Formation Models? Perhaps the 
most serious complication with flow-driven cloud formation 
models is that by themselves they address only one of the 
fundamental timescales of star forming clouds introduced in 
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FIG. 4 Distribution of free-fall times against time in a molecular 
cloud forming in large-scale neutral hydrogen streams. Star forma- 
tion sets in at ~ 11 Myr, when the local free-fall times get substan- 
tially shorter than those in the bulk of the cloud. For comparison, 
substantial CO is visible at ~ 10 Myr. (Adapted from \Heitsch and\ 
\Hartmann\\2008\ .) 



Sec. II.B.l the lag time, tiag, between when clouds begin to 



accumulate and when star formation begins. Taken at face 
value, they do not explain cloud lifetimes, tufo, nor the low 
overall star formation rates or equivalently the long gas de- 
pletion timescales, idcp- This is because these models by 
themselves lack an exit strategy. In the absence of energy 
sources within the cloud, the accumulated mass will start to 
collapse globally, the clouds would settle and convert a sub- 
stantial fraction of their mass into stars ( jHeitsch et a/.||2008b[ 
|Vazquez-Semadeni et al.\ |2007) l, violating not only the ob- 
served cloud lifetime limits, but also the observed limits on 
the star formation rate. Additional processes, most likely stel- 
lar feedback, are required to set the two remaining timescales. 
We come back to this issue in Section ITV] 



C. Modeling Molecular Cloud Formation in 
Hydrodynamic Simulations with Time-Dependent 
Chemistry 

The chemical composition of the ISM is complex. Over 
120 different molecular species have been detected in inter- 



stellar space (Langer et al. 2000 1 and while many of these are 
found in detectable amounts only in dense, well-shielded gas, 
there remain a significant number that have been detected in 
diffuse, unshielded gas ( |0'Neill et a/.||2002| . A full chemical 
model of the ISM can easily involve several hundred different 
atomic and molecular species and isotopologues and several 
thousand different reactions, even if reactions on grain sur- 
faces are neglected, see e.g. the UMIST database ( |Le Teuff] 
|eFar]|2000| ). 

It is currently impractical to incorporate this amount of 
chemistry into a 3D hydrodynamical code. The key to 
constructing time-dependent chemical networks that can run 
alongside the dynamic evolution of the system therefore is 



to select a number of chemical species and mutual reaction 
rates that is small enough so that the chemical network can 
be solved in a short enough time so that it is tractable to do 
so during each system timestep and that is large and complete 
enough so that the overall evolution of the system is still de- 
scribed adequately. In the context of molecular cloud forma- 
tion it is clearly necessary to be able to follow the formation 
and destruction of H2 with a reasonable degree of accuracy. 
Beyond this, the only chemistry that is really required is that 
which plays a role in determining the thermal balance of the 
gas. In other words, we need only follow the chemistry of H2, 
and of a few other major coolants such as C+ or O in low col- 
umn density gas, or CO in high column density regions. As 
few as thirty species and two hundred reactions appear to be 
sufficient for accurately modelling the most important hydro- 
gen, carbon and oxygen chemistry in molecular cloud forma- 
tion calculations ( |Glover et al.\ |2009| l, and a network of this 
size has been shown to be practical to incorporate into a 3D 
hydrodynamical code ( Jappsen et al. 2007| l. 

Many reaction rates are sensitive to the external radiation 
field. Molecular hydrogen, for example, can only remain sta- 
ble in regions where the column density is high enough to 
significantly attenuate the Galactic radiation field fDraine and 
|Bertoldi, 1996). This means, that any sensible calculation of 
chemical reaction rates requires knowledge of the local radia- 
tion field. Ideally, simulations with time-dependent chemistry 
running alongside the hydrodynamics should also include full 
radiative transfer (as discussed in Section [rV.B[ ). This is, how- 
ever, beyond the capabilities of current numerical schemes, 
and most astrophysical applications treat radiation in a very 
approximate fashion only, for example, by assuming a con- 
stant background field or by computing column densities and 
optical depths only along the principle axes of the system. 

Although it has as yet received only limited attention 
from computational astrophysicists, efficient coupling be- 
tween chemical reaction networks and hydrodynamic solvers 
is an active area of research in fields such as combustion mod- 
elling (e.g. Ren and Popej |2007) l or atmospheric chemistry 
(e.g. |Sportisse 2007 1. The basic principles are straightfor- 
ward. One usually uses some form of operator splitting to 
separate the treatment of the chemistry from the advection 
and/or diffusion terms. During the chemistry sub-step, one 
updates the chemical abundances by solving a coupled set of 
rate equations of the form 



~dt 



Ci DiTli 



(3) 



where rii is the number density of species i, and Ci and Di 
are chemical creation and destruction terms that generally de- 
pend on the temperature T and the chemical abundances of 
the other reactants in the system. Most chemical reaction net- 
works are stiff - that is, they contain a wide range of differ- 
ent characteristic timescales - and so to ensure stability, it is 
usually necessary to solve these coupled rate equations with 
an implicit scheme. The simplest implicit techniques have a 
computational cost that scales as the cube of the number of 
species, and so considerable ingenuity has been expended in 
attempts to reduce this cost, for instance by making use of 




FIG. 5 Synthetic emission maps from turbulent-box calculations with time-dependent chemistry. The top left panel depicts the total H2 column 
density, while the top right panel shows the integrated optically thin line emission from the tracer molecule CO. For comparison, the lower 
left image shows the total gas column density. Inspection of the top two images illustrates that CO traces the H2 distribution very well in high 
column density regions, however, fails to do so for the low-density surface layers that are more strongly exposed to the external radiation field. 
This is quantified in the bottom right image which shows the ratio between both values. This ratio is related to the so-called "x-factor" that is 
commonly used to convert CO intensity maps into H2 maps ^Pineda et a/., ,2008^ . (Image courtesy of C. Federrath) 



chemical conservation laws to reduce the number of species 
that need to be tracked, or by taking advantage of the typically 
sparse nature of the Jacobian of the coupled set of equations 
(e.g. |Nejad||2005] l. 



The thermal evolution of the gas is usually modelled us- 
ing a library of cooling functions for each considered species. 
For example, state-of-the-art calculations include the effects 
of atomic fine structure cooling (e.g. C, C"*", O, Si and Si"*"), 
rotational and vibrational cooling of the considered molecules 
(e.g. H2, HD, CO and H2O), Lyman-a cooling, Compton 
cooling, and H+ recombination cooling, as well as other 
processes of lesser importance. The numerical implemen- 
tation usually adopts the isochoric approximation (Springel' 
(2001 J and computes emission strength using the large- 



et al 



velocity gradient approach, which assumes that the emitted 
lines are absorbed locally only. During a given hydrodynamic 
timestep one computes first itad, the rate of change of the in- 



ternal energy due to adiabatic gas physics. Then one has to 
solve an implicit equation for the new internal energy of the 
form 



y„OW ^ yOld 



Wad 



A(p" 



At, 



(4) 



where and u°^^ are the internal energy per unit mass at 
the current and old time, respectively, p""^ is the gas density 
at the current time. It is often necessary to solve this implicit 
equation simultaneously with the chemical rate equations. 
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III. FROM CLOUD CORES TO STARS 

A. Observational Properties of Molecular Cloud Cores 

1 . Statistical Properties 

Emission line observations and dust extinction maps of 
molecular clouds reveal extremely complex morphological 
structure with clumps and filaments on all scales accessible 
by present day telescopes. Typical parameters of different re- 
gions in molecular clouds are listed in Table [l] The volume 
filling factor nj {n) (where n is the local density, while {n) is 
the average density of the cloud) of dense clumps, even denser 
subclumps and so on, is rather small, rangin g from 10% down 



above ^0.5 M0 and a = 1.5 below. The first large study of 



to 0.1% at densities of n > 10^ cm~-^ (|Beuther ef a/.| |2000 



[BliEl [T993i IMcKeel [T9991 [Wilhams ef a/.| |2000| l. This hi 



erarchical configuration is often interpreted as being fractal 
( Bensch et a/.[|200T]|Elmegreen and Falgarone, 1996;,Stutzki, 



et al. 1998|l whi ch however, is still subject to debate ( Blitz and 



Williams T997JI. It is important to note that star formation al- 



ways occurs in the densest regions within a cloud, so only a 
small fraction of molecular cloud matter is actually involved 
in building up stars, while the bulk of the material remains at 
lower densities. This is most likely the key to explaining the 
low star formation efficiencies as discussed above in Section 

urm 

The mass spectrum of clumps in molecular clouds appears 
to be well described by a power law. 



dN 
dm 



(5) 



with the exponent being in the range —1.3 < a < —1.8, 



indicating self-similarity (Kramer et al. 1998 Stutzki and 



|Guesten| |1990| [Williams et al.\ |1994| ). There is no natural 
mass or size scale between the lower and upper limits of the 
observations. The smallest observed structures are protostel- 
lar cores with masses of a few solar masses or less and sizes 
of ^ 0.1 pc. The fact that all studies obtain a similar power 
law is remarkable, and may be the result of turbulent motions 
and thermal instability acting on self-gravitating gas. Given 
the uncertainties in determining the slope, it appears reason- 
able to conclude that there is a universal mass spectrum for the 
clumps within a molecular cloud, and that the distribution is a 
power law within a mass range of three orders of magnitude, 
i.e. from 1 M© to about 1000 M0. Hence, it appears plausible 
that the physical processes that determine the distribution of 
clump masses are rather similar from cloud to cloud. And vice 
versa, clouds that show significant deviation from this univer- 
sal distribution most likely had different dynamical histories. 
Most of the objects that enter in the above morphologi- 



cal analyses are not gravitationally bound (Dib et al. 2007 



Klessen et ar]|2005t|Morata et a/.||20051|Stutzki and Guesten 



1990| l. It is interesting to note that the distribution changes 
as one probes smaller and smaller scales and more and more 
bound objects. When considering prestellar cores, which are 
thought to be the direct progenitors of individual stars or small 
multiple systems, then the mass function is well described by 
a double power law fit dN/dm ex m^" following a — 2.5 



this kind was published by Motte et al. ( 1998), for a popula- 
tion of submillimetre cores in p Oph. Using data obtained 
with IRAM, they discovered a total of 58 starless clumps, 
ranging in mass from 0.05 Mq to 3 Mq. Similar re- 
sults are obtained from the Serpens cloud (Testi and Sargentj 
|1998| ), for Orion B North Pohnstone et al., ,2001. ) and Orion 
B South (Johnstone et al. 2006| l, or for the Pipe Nebula (Lada| 
'et all '2006'). Currently all observational data ( Alves et al. 
2007; Di Francesco et a/.!l2007!lJohnstone et a/.l 12001 2006 
2000; Lada et al. 2008 , Motte et al. 1998 ; Nutter and Wari_ 
Thompson 2006; |Testi and Sargent^ |1998f [Ward-Thompson 



.et al , 2007 ) reveal that the mass function of prestellar cores 
is strikingly similar in shape to the stellar initial mass func- 
tion, the IMF. To reach complete overlap one is required to 
introduce a mass scaling or efficiency factor in the range 2 
to 10, which differs in different regions. An exciting inter- 
pretation of these observations is that we are witnessing the 
direct formation of the IMF via fragmentation of the parent 
cloud. However, we note that the observational data also in- 
dicate that a considerable fraction of the prestellar cores do 
not exceed the critical mass for gravitational collapse, much 
like the clumps on larger scale. The evidence for a one-to-one 
mapping between prestellar cores and the stellar mass thus is 
by no means conclusive as we will discuss in more detail in 
Section UlLCl 



2. Individual Cores 

Density Structure. The density structure of prestellar cores 
is typically estimated through the analysis of dust emission 
or absorption using near-IR extinction mapping of back- 
ground starlight, mapping of millimeter/submillimeter dust 
continuum emission, and mapping of dust absorption against 
the bright mid-IR background emission (Bergin and Tafalla] 



2007 ). A main characteristic of the density profiles derived 
with the above techniques is that they require a central flat- 
tening. The density gradient of a core is flatter than r~^ 
within radii smaller than 2500 - 5000 AU, and th at the typ- 



ical central density of a core is 10 - 10 cm ( .Motte et al. 



1998 Ward- Thompson et al. 1999| ). A popular approach is to 
describe these cores as truncated isothermal (Bonnor-Ebert) 
sphere ( |Bonnor| |1956[ |Ebert| |1955| ), that often (but not al 



ways) provides a good fit to the data ( Alves et al. 2001 Bac 



mann et al. 2001 Kandori et al. 2005 ). These are equilib- 
rium solutions of self-gravitating gas spheres bounded by ex- 
ternal pressure. However, such density structure is not unique. 
Numerical calculations of the dynamical evolution of super- 
sonically turbulent clouds show that transient cores forming at 
the stagnation points of convergent flows exhibit similar mor- 



phology (Ballesteros-Paredes et al. 2003 1. 

Thermal Stucture. The kinetic temperature of the dust 
and gas components in a core is regulated by interplay be- 
tween various heating and cooling processes. At high densi- 
ties (> 10'' cm^'^) in the inner part of the cores, the gas and 



dust have to be coupled thermally vi a collisions ([Burke and| 
HollenbachI [1983] [GoidsrmQil [Goldsmith and Langer| 
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t = 22-50 Myr 




FIG. 6 Formation and growth of molecular cloud cores by thermal instability triggered by a large-scale convergent flow: A small cold 
condensate grows from the thermally unstable warm neutral medium by outward propagation of its boundary layer Coalescence and merging 
with nearby clumps further increases its mass and size. The global gravitational potential of the proto-cloud enhances the merging probability 
with time. The images show 2D slices of the density (logarithmic colour scale) and the gas velocity (indicated as arrows) in the plane 
perpendicular to the large scale flow. (From a numerical simulation by \Banerjee et aL| ( [200Sp ) 



|1978V At lower densities, which correspond to the outer parts 
of the cores, the two temperatures are not necessary expected 
to be the same. Thus, the dust and gas temperature distri- 
butions need to be inferred from observations independently. 
Large-scale studies of dust temperature show that the grains in 
starless cores are colder than in the surrounding lower-density 
medium. Far-IR observations toward the vicinity of a number 
of dense cores provide evidence for flat or decreasing tem- 
perature gradients with cloud temperatures of 15 — 20 K and 
core values of 8 - 12 K ( jToth et al.\ [2004} |Ward-Thompson| 
\et al.\ |2002[ ). These observations are consistent with dust ra- 
diative transfer modeling in cores illuminated by interstellar 
radiation field (|Keto and FieIdl|2005||Langer et a/.||2005||StF 



|matellos et al.\ 2007b| l, where the dust temperature is ~ 7 K 
in the core center and increases up to 16K in the envelope. 
The gas temperature in molecular clouds and cores is com- 
monly infered from the level excitation of simple molecules 
likeCO and NH3 ( |Evans| [7999} [Walmsley and Ungerechts} 



1983 I. One finds gas temperatures of 10-15 K, with a pos- 



sible increase toward the lower density gas near the cloud 
edges. It is believed that the gas heating in prestellar cores 
mostly occurs through ionization by cosmic rays, while the 
cooling is mainly due to line radiation from molecules, es- 
pecially CO (Goldsmith and Langer 1978| l. Altogether, the 
fact that prestellar cores are cold and roughly isothermal with 
at most a modest increase in temperature from the center to 
the edge is consistent with numerical models of cores form- 
ing from thermal instab ility (|Banerjee| [2008 [ [Heitsch et al.\ 
|2006a| |Keto and CaselEl |2008p , see also Figure |6] 

Chemical Stucture. Maps of integral line intensity can look 
very different for different molecular tracers. In particular, the 
N2H+ and NH3 emission more closely follows the dust emis- 
sion while the C^*0 and CS emission appears as a "ring-like" 
structure around the dust emission maximum ([Bergin et at. 



2OO2I ILada et al.\ |2003i |Maret et al.\ |2007l [Tafalla et at. 
2002|l. For illustration see Figure [7] The common theoreti 



cal interpretion of these data is that carbon-bearing species, 
represented by CO and CS freeze-out on the dust grains from 
the gas while the abundances of nitrogen-hydrogen bearing 



molecules, N2H"'" and NH3, either stay constant or decay 
more slowly. At the same time, chemical models of prestel- 
lar cores predict that molecules in the core envelope have to 



be destroyed by interstellar UV field (Aikawa et at. 2008 



[Pavlyuchenkov et al. 2006 1. The chemical stratification sig 
nificantly complicates the interpretation of molecular fine ob- 
servations and again requires the use of sophisticated chemical 
models which have to be coupled to the dynamical evolution. 
From observational side, the freeze-out of many molecules 
makes it difficult to use their emission lines for probing the 
physical conditions in the inner regions of the cores. At the 
same time, the modeling of the chemical evolution can pro- 
vide us with the important parameters of the cores. For exam- 
ple, the level of CS depletion can be used to constrain the age 
of the prestellar cores while the deficit of CS in the envelope 



can indicate the strength of the external UV field ( jBergin and 
|Tafalla[[2007| l. In any case, any physical interpretation of the 
molecular lines in prestellar cores has to be based on chemi- 
cal models and should do justice to the underlying density and 
velocity pattern of the gas. 

Kinematic Stucture. In contrast to the supersonic velocity 
fields observed in molecular clouds, dense cores have low in- 
ternal velocities. Starless cores in clouds like Taurus, Perseus, 
and Ophiuchus systematically present spectra with close-to- 
thermal linewidths, even when observed at low angular reso- 



lution pijina et a/.]|1999[|Myers||1983| l. This indicates that the 
gas motions inside the cores a re subsonic or at best transsonic, 
i.e. with Mach numbers ^ 2 (Andre et al. 2001 ; Kirk et al. 



|2007t [Rosolowsky et al.\ |2008| l. In some cores also inwarc 
motions have been detected. They are inferred from the ob- 
servation of optically thick, self-absorbed lines of species like 
CS, H2CO, or HCO+, in which low-excitation foreground gas 
absorbs part of the background emission. Typical inflow ve- 
locities are of order of 0.05 — O.lkm/s and are observed on 
scales of 0.05 — 0.15pc, comparable to the observed size of 



the cores (Lee et al. 1999j ). The overall velocity structure 
of starless cores appears broadly consistent with the structure 
predicted by models in which protostellar cores form as the 
stagnation points of convergent flows, but the agreement is not 
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FIG. 7 Maps of molecular line emission from C^*0, N2H+, and CS superimposed on a dust extinction maps of the dark cloud core Barnard 
68 jAlves ef a/.[[2001 , Bergin et al. , 2002 , Lada et a/.[|2003^ . The three images illustrate the effects of depletion onto grains in the high-density 
central region of the core. N2H"'^ is the least and CS the most depleted species. (Image courtesy ofE. A. Bergin) 



perfect. Simulations of core formation do correctly find that 
most cores are at most transsonic ( |Klessen et a/.| |2005 ; Offner] 



et al. 2008 ), but the distribution of velocity dispersions has 



a small tail of highly supersonic cores that is not observed. 
Clearly more theoretical and numerical work is needed. In 
particular, the comparison should be based on synthetic line 
emission maps, which requires to couple a chemical network 
and radiative transfer to the simulated density profiles as dis- 
cussed above. In addition, it is also plausible that the discrep- 
ancy occurs because the simulations do not include all the nec- 
essary physics such as radiative feedback and magnetic fields. 
Subsonic turbulence contiibutes less to the energy budget of 
the cloud than thermal pressure and so cannot provide suffi- 



cient support against gravitational collapse ( [Goodman et al. 



[1998; Myers, 1983; Tafal la et a/.||2006| . If cores are longer 
lasting entities there must be other mechanisms to provide sta- 
bility. Obvious candidates are magnetic fields (Shu et al. 
1987| l. However, they are usually not strong enough to pro 



vide sufficient support ( Bourke et al.' 20011 Crutcher et al. 
|2008bt [Crutcher and Troland 2000, Crutcher et al.\ |1999 
as discussed below. Most observed cores are thus likely to 
be evolving transient objects that never reach any equilibrium 
state. 

Magnetic Field Structure. Magnetic fields are ubiquitously 
observed in the interstellar gas on all scales ( jCrutcher et al.\ 
|2003[ [Heiles and Troland] |2005[ ). However, their importance 
for star formation and for the morphology and evolution of 
molecular cloud cores remains controversial. A crucial pa- 
rameter in this debate is the ratio between core mass and mag- 
netic flux. In supercritical cores, this ratio exceeds a criti- 
cal value and collapse can proceed. In subcritical ones, mag- 



gin but very often by factors of many if non-detections are 
include d ( |Bourke et a/1|200T)[CruFcher||1999[|Crutcher et al.\ 
2008b[ l. The polarization of dust emission offers an alternative 



pathway to studying the magnetic field structure of molecu 
lar cloud cores. MHD simulations of turbulent clouds pre- 
dict degrees of polarization between 1 and 10%, regardless of 
whether turbulent energy dominates over the magnetic energy 
(i.e. the turbulence is super- Alfvenic) or not ( |Padoan et al.\ 
|2001[[Padoan and Nordlund||1999| l. However, converting po- 
larization into magnetic field strength is very difficult ( Heitsch 
et al. 200 lb). Altogether, the current observational finding 



imply that magnetic fields must be considered when studying 
stellar birth, but also that they are not the dominant agent that 
determines when and where stars form within a cloud. Mag- 
netic fields appear too weak to prevent gravitational collapse 
to occur 

This conclusion means that in many cases and to rea- 
sonable approximation purely hydrodynamic calculations are 
sufficient for star formation simulations. However, when 
more precise and quantitative predictions are desired, e.g. 
when attempting to predict star formation timescales or bi- 
nary properties, it is necessary to perform magnetohydrody- 
namic (MHD) simulations or even consider non-ideal MHD. 
The latter means to take ambipolar diffusion (drift between 
charged and neutral particles) or Ohmic dissipation into ac- 
count. Recent numeiical simulations have shown that even a 
weak magnetic field can have noticeable dynamical effects. It 
can alter how cores fragment ( Hennebelle and Fromang 200^ 
[Hennebelle and Teyssief] [2008,, Price and Bate„2007b,.2008p , 
change the coupling between stellar feedback processes and 



their parent clouds (Krumholz et al. 2007c Nakamura and 



netic fields provide stability (Mouschovias 1991a b ; ,Spitzer Li 2007j), influence the properties of protostellar disks due to 



|1978| l. Measurements of the Zeeman splitting of molecular 
lines in nearby cloud cores indicate mass-to-flux ratios that lie 
above the critical value, in some cases only by a small mar- 



magnetic braking ( [Mellon and Li| [2008a|bt [Price and Bate 



2007a I, or slow down the overall evolution (Heitsch et al. 



2001a I 
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3. Models of Cloud Evolution and Star Formation 

There are two main competing models that describe the 
evolution of the cloud cores. It was proposed in the 1980's 
that cores in low-mass star-forming regions evolve quasi- 
statically in magnetically subcritical clouds ( |Shu ef ar||1987| l. 
Gravitational contraction is mediated by ambipolar diffusion 
(|Mouschovias| 1 1 976| [T97 9 , 1991ai Mous chovias and Paleolo-| 
|gou |1981| l causing a redistribution of magnetic flux until the 
inner regions of the core become supercritical and go into dy- 
namical collapse. This process was originally thought to be 
slow, because in highly subcritical clouds the ambipolar diffu- 
sion timescale is about 10 times larger than the dynamical one. 
However for cores close to the critical value, as is suggested 
by observations, both timescales are comparable. Numerical 
simulations furthermore indicate that the ambipolar diffusion 
timescale becomes significantly shorter for turbulent veloci- 
ties similar to the values observed in nearby star-forming re- 



gion (jFatuzzo and Adams |2002[ peitsch et al\ |2004[ | Li and| 



[Nakamura |2004| l. The fact that ambipolar diffusion may not 
be a slow process under realistic cloud conditions, as well 
as the fact that most cloud cores are magnetically supercrit- 



ical (Bourke et al. 


2001; Crutcher et al.\ |2008bt |Crutcher| 


land Troland, ,2000, 


Crutcher et a/. [11999 1 has cast significant 



doubts on any magnetically-dominated quasi-static models of 
stellar birth. 

For this reason, star-formation research has turned into con- 
sidering supersonic turbulence as being on of the principal 
physical agent regulating stellar birth. The presence of tur- 
bulence, in particular of supersonic turbulence, has impor- 
tant consequences for molecular cloud evolution. On large 
scales it can support clouds against contraction, while on 
small scales it can provoke localized collapse. Turbulence 
establishes a complex network of interacting shocks, where 
dense cores form at the stagnation points of convergent flows. 
The density can be large enough for gravitational collapse to 
set in. However, the fluctuations in turbulent velocity fields 
are highly transient. The random flow that creates local den- 
sity enhancements can disperse them again. For local col- 
lapse to actually result in the formation of stars, high den- 
sity fluctuations must collapse on timescales shorter than the 
typical time interval between two successive shock passages. 
Only then are they able to 'decouple' from the ambient flow 
and survive subsequent shock interactions. The shorter the 
time between shock passages, the less likely these fluctua- 
tions are to survive. Hence, the timescale and efficiency of 
protostellar core formation depend strongly on the wavelength 
and strength of the driving source (Ballesteros-Paredes et al.\ 
20071 'Heitsch eTaT] [2001a) [Klessen et al. . 2000 ; .Krumholz 



and McKee 2005; Mac Low and Klessen, 2004; McKee and 



Ostriker[ |2007[ [Vazquez-Sema deni et al. 2003 ), and accre- 



tion histories of individual protostars are strongly time vary- 



ing (Klessen 2001b Schmeja and Klessen 2004| i. 

Interstellar turbulence is observed to be dominated by large- 
scale modes ( Mac Low and Ossenkopf[ 2000, pssenkopf[ 
\et al.\ [200T] [Ossenkopf and Mac Lowj [2002[ l. This implies 
it is very efficient in sweeping up molecular cloud material, 
thus creating massive coherent structures. The result is a large 



region in which many Jeans masses of material become un- 
stable to collapse at about the same time, leading to coherent 
structures in the forming stars. This is a likely explanation for 



the observed clustering of young stars (Klessen et al. 20001, 
as we discuss in the following section. 



B. Spatial Distribution 

The advent of sensitive infrared detectors in the last decade 
has made it possible to perform wide-area surveys. These 
have led us to recognize that most stars form in clusters and 
aggregates of various size and mass scales, and that isolated 
or widely distributed star formation is the exception rather 
than the rule (Lada and Lada 2003 [ l. The complex hierar- 
chical structure of molecular clouds (Figure |2]) provides a 
natural explanation for this finding. Star-forming molecular 
cloud cores vary enormously in size and mass. In small, low- 
density, clouds stars form with low efficiency, more or less in 
isolation or scattered around in small groups of up to a few 
dozen members. Denser and more massive clouds may build 
up stars in associations and clusters of a few hundred mem- 
bers. This appears to be the most common mode of star for- 
mation in the solar neighborhood ( Adam s and Myers] [2001 [ I . 
Examples of star formation in small groups and associations 



are found in the Taurus-Aurigae molecular cloud (Hartmann 
[2002[ ). Young stellar groups with a few hundred members 
form in the Chamaeleon I ([Persi et al.\ |2000| or p-Ophiuchi 
(Bontemps 2001) dark clouds. Each of these clouds is at a 



distance of about 130 to 160pc from the Sun. Like most of 
the nearby young star forming regions they appear to be as- 
sociated with a ring-like structure in the Galactic disk called 



Gould's belt ( [Poppel[[T997 1. 

The formation of dense rich clusters with thousands of 
stars is rare. The closest region where this happens is the 
Orion Nebula Cluster in L 1641 ( |Hillenbrand| [T^ [Hiller> 



brand and Hartmann 1998 1, which lies at a distance of 410 pc 
(Caballero 2008'; 'Hir ota eFall [2007) [Menten et aL\ [2057j 
iSandstrom et a/., ,2007) . A rich cluster somewhat further away 



is associated with the Monoceros R2 cloud ( Carpenter et al. 
1997| ) at a distance of - 830 pc. The cluster NGC 3603 is 
roughly ten times more massive than the Orion Nebula Clus- 
ter It lies in the Carina region, at about 7kpc distance. It 
contains about a dozen O stars, and is the nearest object anal- 



ogous to a starburst knot (Brandl et al. 1999 Moffat et al. 



[2002[ ). To find star-forming regions building up hundreds of O 
stars one has to look towards giant extragalactic Hll regions, 
the nearest of which is 30 Doradus in the Large Magellanic 
Cloud, a satellite galaxy of our Milky Way at a distance at 
55 kpc. The giant star forming region 30 Doradus is thought 
to contain up to a hundred thousand young stars, including 



more than 400 O stars (Hunter et al. 1995 Townsley et al. 



[2006; ,Walborn and Blades, ,1997 1. This sequence as depicted 
in Figure[8[demonstrates that the star formation process spans 
many orders of magnitude in scale, ranging from isolated sin- 
gle stars to massive young clusters with several 10^ stars. 
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FIG. 8 Comparison of clusters of different masses scaled to same 
relative distance. The cluster in the upper left comer is the Orion 
Nebula Cluster (IMcCaughrean 200 l| l and the one at the lower left is 
NGC 3603 (Bran dl et all 1 1999 ; both observed with the Very Large 



Telescope at infrared wavelength. The large cluster in the center is 
30 Doradus in the LMC observed with the Hubble Space Telescope 
(courtesy of M. J. McCaughrean). The total mass increases roughly 
by a factor of ten from one cluster to the other. { Composite image 
courtesy ofH. Zinnecker) 



C. The Stellar initial Mass Function and other Statistical 
Characteristics of Star Formation 

The mass distribution of young stars follows a well-known 
distribution called the Initial Mass Function (IMF). For stellar 
masses m > IMq it shows a power-law behavior dN/dm cx 
m" , with slope a = -2.3 ( |Chabrier| 2003; Kr oupa) [2002{ 



Salpeterj 11955] [Scalol [T998] l. Below IA/q, the IMF flat- 



tens, a change in behavior than can be represented either as 
a lognormal (Chabrier 2003| l or a change in power law in- 



dex (Kroupa 2002 1. At the extreme ends of the stellar mass 
spectrum, however, our knowledge of both the IMF are lim- 
ited. Massive stars are very rare and rather short lived. The 
number of massive stars that are sufficiently near to study in 
detail and with very high spatial resolution, for example to de- 
termine multiplicity, therefore is small ( [Zinnecker and Yorke[ 
|2007| l. Low-mass stars and brown dwarfs, on the other hand, 
are faint, so they too are difficult to study in detail ( [Burrows 



et al. 2001 1. Such studies, however, are in great demand, 
because secondary indicators such as the fraction of binaries 
and higher-order multiples as function of mass, or the distri- 
bution disks around very young stars or possible signatures of 
accretion during their formation are probably better suited to 
distinguish between different star-formation models than just 
looking at the IMF. In contrast to the observational agreement 
on the IMF, at least above the substellar regime, there is still 



considerable disagreement on the theoretical side. The ori- 
gin of the IMF is a major topic of theoretical research which 
we examine only briefly here to give the necessary theoretical 
background for our discussion of numerical work. Other re- 



views provide considerably more detail (Bonne II et al. 2007 
Larson 2007J |Mac Low and KIessen| |2004[ |McKee and Os-| 
5rker,.2007| . 



Early models for the origin of the IMF generally relied on 
statistical arguments, appealing to random process of collapse 
in a fractal cloud ( ,Elmegreen[ |1997[ |2000| l, or to the cen- 
tral limit theorem to explain its characteristic shape ( |Larson[ 
|I973[ [Zinnecker! [1984 '). Researchers have also invoked feed- 
back processes that cut off accretion onto individual protostars 
( [Adams and Fatuzzo||1996| . Today, however, there are three 
dominant schools of thought regarding the origin on the IMF, 
although the boundaries between these pictures are not clearly 
defined, and numerous hybrid models have been proposed. 

One model, called core accretion, takes as its starting point 
the striking similarity between the shape of the observed core 
mass distribution and the IMF. This model posits that there 
is a one-to-one relation between the distributions, so that in- 
dividual cores are the progenitors of individual stars or star 
systems. The factor of ^ 3 decrease in mass between cores 
and stars is the result of feedback processes, mostly protostel- 
lar outflows, that eject a fixed fraction of the mass in a core 
rather than letting it accrete onto the star ( Matzner and Mc-] 
Kee 2000| l. This model reduces the problem of the origin of 



the IMF to the problem of determining the mass spectrum of 
bound cores, although strictly speaking the idea that the IMF 
is set by the mass spectrum of cores is independent of any par- 
ticular model for the origin of that mass spectrum. Arguments 
to explain the core mass distribution generally rely on the 
statistical properties of turbulence (H ennebelle an d Chabrierj 
|2008; Klessen, 2001ai [Padoan and N ordlund 2002), which 
generate structures with a pure powerlaw mass spectrum. The 
thermal Jeans mass in the cloud then imposes the flattening 
and turn-down in the observed mass spectrum. 

A second model for the origin of the IMF, called compet- 
itive accretion, focuses instead in interaction between proto- 
stars, and between a protostellar population and the gas cloud 



around it (jBate and BonneU 2005||Bate et al. 


2003 




BonneIl| 


[and Bate[ j2002| Bonnell et al. 2001a [2008 


2006 




2001b 1. 



In the competitive accretion picture the origin of the peak in 
the IMF is much the same as in the core accretion model: it 
is set by the Jeans mass in the prestellar gas cloud. How- 
ever, rather than fragmentation in the gas phase producing a 
spectrum of core masses, each of which collapses down to a 
single star or star system, in the competitive accretion model 
all gas fragments down to roughly the Jeans mass. Prompt 
fragmentation therefore creates a mass function that lacks the 
powerlaw tail at high masses that we observe in the stellar 
mass function. This part of the distribution forms via a sec- 
ond phase in which Jeans mass-protostars compete for gas in 
the center of a dense cluster The cluster potential channels 
mass toward the center, so stars that remain in the center grow 
to large masses, while those that are ejected from the clus- 



ter center by iV-body interactions remain low mass (Bonnell 



et a/.[[2004[ [Klessen and Burkert(,2000j . In this model, the ap- 
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parent similarity between the core and stellar mass functions 
is an illusion, because the observed cores do not correspond 
to gravitationally bound structures that will collapse to stars 
HUark and Boiinelll [20061 [Smith et a/.||2008] l. 

One potential drawback to both the core accretion and com- 
petitive accretion models is that they rely on the Jeans mass 
to determine the peak of the IMF, but leave unanswered the 
question of how to compute it. This question is subtle be- 
cause molecular clouds are nearly isothermal, but they con- 
tain a very wide range of densities, and it is unclear which 
density should be used. A promising idea to resolve this ques- 
tion, which is the basis for a third model of the IMF, focuses 
on the thermodynamic properties of the gas. The amount of 
fragmentation occurring during gravitational collapse depends 



on the compressibility of the gas ( Li et al. 2003 1. For poly- 
tropic indices 7 < 1, turbulent compressions cause large den- 
sity enhancements in which the Jeans mass falls substantially, 
allowing many fragments to collapse. Only a few massive 
fragments get compressed strongly enough to collapse in less 
compressible gas though. In real molecular gas, the compress- 
ibility varies as the opacity and radiative heating increase. 
Larson ( 2005| l noted that the thermal coupling of the gas to 



the dust at densities above riciit ^ 10 — 10 cm^'' leads to 
a shift from an adiabatic index of 7 ~ 0.7 to 1.1 as the den- 
sity increases above ricrit- The Jeans mass evaluated at the 
temperature and density where this shift occurs sets a mass 
scale for the peak of the IMF. The apparent universality of the 
IMF in the Milky Way and nearby galaxies may be caused by 
the insensitivity of the dust temperature on the intensity of the 
interstellar radiation field (Elmegre en et al.\ 2008[ ). Not only 
does this mechanism set the peak mass, but also appears to 
produce a power-law distribution of masses at the high-mass 
end comparable to the observed distribution ( [Jappsen et al.\ 
|2005l l. 

Each of these models has potential problems. In the core 
accretion picture, hydrodynamic simulations seem to indicate 
that massive cores should fragment into many stars rather than 



collapsing mon olithically (,Bonnell and Bate, 2006 Clark and 



|Bonnell| |2006[ |Dobbs et al. 2005| l. The hydrodynamic sim 
ulations almost certainly suffer from over-fragmentation be- 
cause they do not include radiative feedback from embed- 
ded massive stars ( Krumholz 2006j |Krumholz et al.\ |2007a[ 
IKrumholz and McKee, |2008j), but no simulation to date has 
successfully formed a massive core in a turbulent cloud and 
followed it all the way to the formation of a massive star 

In addition, the suggestion of a one-to-one mapping be- 
tween the observed clumps and the final IMF is subject to 
strong debate. Many of the prestellar cores discussed in Sec- 
appear to be stable entities ( [Johnstone et a/.|[2001[ 



tion 



III.A.1 



|2006||2000 



[Lada et aT] |2008| l, and thus are unlikely to be in 
a state of active star formation. In addition, the simple in- 
terpretation that one core forms on average one star, and that 
all cores contain the same number of thermal Jeans masses, 
leads to a timescale problem d Clark et oT] |2007| ) that requires 
differences in the core mass function and the IMF. 

The criticism regarding neglect of radiative feedback ef- 
fects also applies to the gas thermodynamic idea: the cooling 
curves which that model assumes in order to derive the Jeans 



mass ignore the influence of protostellar radiation on the tem- 
perature of the gas, which simulations show can suppress frag- 
mentation in at least some circumstances ( jKrumholz et al.\ 



2007a I. The competitive accretion picture has also been chal- 



lenged, on the grounds that the kinematic structure observed 
in star-forming regions is inconsistent with the idea that pro- 
tostars have time to interact with one another strongly before 



they completely accrete their parent cores (Andre et al. 2007 
[Krumholz eFaT] |2005c1 ). 



D. Modeling Cloud Fragmentation and Protostellar 
Collapse 

To adequately model the fragmentation of molecular 
clouds, the formation of dense cloud cores, the collapse of the 
gravitationally unstable subset of cores, and finally the build- 
up and mass growth of embedded protostars in their interior 
is an enormous computational challenge. It requires to fol- 
low the evolution of self-gravitating, highly turbulent gas over 
many order of magnitudes in density and lengthscale. Ow- 
ing to the stochastic nature of supersonic turbulence, it is not 
know in advance where and when local collapse occurs. One 
therefore needs highly flexible numerical methods for solv- 
ing the equations of hydrodynamics, schemes that can pro- 
vide sufficient degrees of precision and resolution throughout 
the entire computational domain in an adaptive fashion. 

The star formation community is following two highly 
complementary approaches to reach these goals. One set of 
methods is based on dividing the computational domain into 
small volume elements and follow the fluxes of all relevant 
quantities from one cell to the other. These grid-based meth- 
ods adopt an Eulerian point of view, because the flow is fol- 
lowed from fixed positions in space. A popular alternative is 
to split the model cloud into individual parcels of gas and fol- 
low their mutual interaction and evolution. Particle-methods 
therefore correspond to a Lagrangian point of view following 
the trajectories of individual fluid elements. 



1 . Grid-Based Methods 

The mathematical formulation of hydrodynamics consists 
of a set of partial differential equations that relate different 
flow properties (such as density and velocity) with each other 
and with thermodynamic quantities (e.g. pressure, tempera- 
ture or internal energy of the medium). They can be formu- 
lated in conservative form corresponding to the conservation 
of mass, momentum, and energy. As the number of variables 
is larger than the number of equations in the system, an ad- 
ditional equation is needed to find a unique solution. This 
closure relation is called equation of state and usually speci- 
fies the pressure as function of other thermodynamic variables 
( [Landau and Lifshitz] [1983[ ). In a broad sense, the hydrody- 
namics equations describe how signals propagate through a 
medium. They specify how local quantities relate to fluxes, 
e.g. how the density in some control volume depends on the 
mass flux through its surface. Equations of this kind are called 
hyperbohc equations. 
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Numerical solutions to partial differential equations always 
require discretization of the problem. This means that instead 
of continuous space and time dimensions we consider a dis- 
crete set of points. The computational domain is subdivided 
into individual volume elements surrounding node points on a 
grid or unstructured mesh. Finite volume methods are proce- 
dures for representing and evaluating partial differential equa- 
tions as algebraic equations. They play a key role in compu- 
tational fluid dynamics. Similar to finite difference schemes, 
values are calculated at discrete places on a meshed geom- 
etry. Volume integrals that contain a divergence term are 
converted to surface integrals, using the divergence theorem. 
These terms are then evaluated as fluxes at the surfaces of each 
volume element. Because the flux entering a given volume 
is identical to that leaving the adjacent volume, these meth- 
ods are conservative. Finite volume methods have been in the 
focus of applied mathematics for decades. They have well 
defined convergence properties and available code packages 
have reached a very high degree of maturity and reliability. 

Finite volume schemes are most easily formulated for 
Cartesian grids with fixed cell size. The cell size determines 
the spatial resolution of the code. Wherever higher resolution 
is needed, it can be achieved by refining the grid. We speak 
of adaptive mesh refinement (AMR), when this is done in an 
automated and locally adjustable way. There are a number of 



different approaches to AMR in the literature (Plewa et al. 
[2005.) . Most AMR treatments are based on finite-element 
models on unstructured meshes. They have the advantage 
to adapt easily to arbitrary complicated boundaries, however, 
constructing the mesh is very time consuming. When using 
Cartesian grids, one can refine on individual cells or on larger 
groups of cells, so-called blocks (Bell et al\ |1994[ |Berger| 



and ColeUa 1989; Berger and Ohger 19841. Cartesian AMR 



codes nowadays belong to the standard repertoire of numeri- 
cal star formation studies. 

In the following we list a few popular hydrodynamic and 
magnetohydrodynamic codes that have been developed in the 
past decade. All but two are freely available, although in 
some cases registration is needed before being able to down- 
load it from the web. ZEUS-MP (http://cosmos.ucsd.edu/lca- 
yyww/ software/index.html) is a parallel, non-adaptive hydro- 
and magnetohydrodynamics code with self-gravity and ra- 
diation ("Hayes et al. 



20061 |Norman| [2000) l. NIR 



VANA (htpp://nirvana-code. aip.de) is an AMR code for non- 
relativistic, compressible, time-dependent, ideal or nonideal 
(viscosity, magnetic diffusion, thermal conduction) MHD 
FLASH (http://flash.uchicago.edu/) is a 



( |Ziegler| [2005]). 
highly modular, parallel adaptive-mesh code initially de- 
signed for thermonuclear runaway problems but also capable 
of a wide variety of astrophysical problems ( [Fryxell et al.\ 
[2000). ENZO (http://lca.ucsd.edu/projects/enzo) is a hy- 
brid AMR code (hydrodynamics and A^-body) which is de- 
signed to do simulations of cosmological structure formation 



(O'Sheae/a/. 2004). It has been extended to include mag- 
netic fields, star formation, and ray-tracing radiation trans- 
fer. ATHENA ( |http : //w w w.astro .princeton. edu/j stone/ athena| l 
is an MHD code built on a flexible framework that is designed 
to allow easy and modular extension to include a wide vari- 



ety of physical processes (IStone et al. 2008[ l. The public ver- 
sion is non-adaptive, contains only hydrodynamics and MHD, 
and uses a fixed Cartesian grid, but extensions exist for non- 
Cartesian grids, for static and adaptive mesh refinement, for 



gravity, and for ionizing radiative transfer (Krumholz et al. 



|2007c| l. Another widely used AMR code for MHD calcula- 
tions is RAMSES ( |Fromang et'al\ [20061 [Teyssier| [2002| . It 
is very versatile with applications ranging from the the two- 
phase interstellar medium, to star and planet formation, as 
well as cosmological structure formation. A grid code use 
commonly in star formation simulations, but which is not 
publicly available, is ORION: a parallel hydrodynamics code 
that includes self-gravity, sink particles coupled to a proto- 
stellar evolution model, and diffuse radiative transfer ( Fisherj 



2002 



et al 



Kieiiil [T9991 [Krumholz et oT] |2007b| [2004j |Truelove| 



1998[ l. Last in our short summary is Proteus, a finite- 



volume method based on a gas-kinetic formulation of the 
microscopic transport properties (|Prendergast and Xul [1993} 



Slyz and Prendergast 1999 Xu 2001 1. This approach al 



lows the user to fully control the dissipative effects, making 
the scheme very attractive for e.g. turbulent transport studies. 
However, adding additional physics is generally more compli- 
cated than for other schemes. Proteus is fully parallelized and 
includes self-gravity, magnetic fields and a two-fluid model 
for ambipolai- drift ( [Heitsch et a/:i[2008bl [20071 [2004l l. 



2. Particle-Based Methods 

Using a particle based scheme to solve the equations of 
hydrodynamics was first introduced by Lucy ( 1977[) an d pro- 
posed independently by [Gingold and Monaghan ( 1977[ ). Orig- 
inally envisioned as a Monte-Carlo approach to calculate the 
time evolution of a hydrodynamic system, the formalism of 
smoothed particle hydrodynamics (SPH) is more intuitively 
understood as a particle interpolation scheme ( [Gingold and| 
Monaghan, 1982). This provides better estimates for the er- 
rors involved and the convergence properties of the method. 
Excellent overviews of the method and some of its appli- 
cations provide the reviews of |Benz[ ( |1990[ ) and [Monaghan] 
( [T992l[2U05] l. 

In the framework of classical physics, fluids and gases are 
large ensembles of interacting particles with the state of the 
system being described by the probability distribution func- 
tion in phase space. Its time evolution is governed by Boltz- 
mann's equation ( [Landau and Lifshitz 1983| l. Hydrodynamic 
quantities can then be obtained in a local averaging process 
involving scales larger than the local mean-free path. A re- 
lated approach is facilitated in SPH. The fluid is represented 
by an ensemble of particles i, each carrying mass, momentum, 
and hydrodynamical properties. The technique can therefore 
be seen as an extension to the well known TV-body methods 
used in stellar dynamics. Besides being characterized by its 
mass iTii and velocity Vi and its location r^, each particle is 
associated with a density pi, an internal energy (equivalent 
to a temperature T^), and a pressure pi. The time evolution 
of the fluid is then represented by the time evolution of the 
SPH particles. Their behavior is governed by the equation of 
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motion, supplemented by further equations to modify the hy- 
drodynamical properties. Thermodynamical observables are 
obtained by averaging over an appropriate subset of the SPH 
particles. 

Mathematically, the local averaging process for a quantity 
/(r) can be performed by convolution with an appropriate 
smoothing function W{r,h): 



(fir)) 



f{r')W{r-r',h)d^r' . 



(6) 



This function W{r, h) is often referred to as the smoothing 
kernel. It must be normalized and approach the Dirac delta 
fuction in the limit h — > 0. For simplicity, most authors 
adopt spherical symmetry in the smoothing and averaging pro- 
cess, i.e. the kernel degrades to an isotropic function of the 
interparticle distances: W{r,h) = W{r,h) with r = |r | and 
\h\. 

The basic concept of SPH is a particle representation of the 
fluid. Hence, the spatial integration in the averaging process 
transforms into a summation over a fixed number of points. 
For example, the density at the position of particle i is com- 
puted as 



\,h) 



(7) 



In this picture, the mass of each particle is smeared out over 
the kernel region. The continuous density distribution of the 
fluid is then obtained by summing over the local contribution 
from neighboring elements j. The name "smoothed particle 
hydrodynamics" derives from this analogy. 

In star formation studies SPH is popular because it is in- 
trinsically Lagrangian. As opposed to mesh-based methods, it 
does not require a fixed grid to represent fluid properties and 
calculate spatial derivatives ( [Hockney and Eastwood[ |1988| l. 
The fluid particles are free to move and - in analogy - con- 
stitute their own grid. The method is therefore able to re- 
solve very high density contrasts, by increasing the particle 
concentration where needed. This it most effective, if the 
smoothing length is adaptable (Monaghan 2005| l. There is 
no need for the complex and time-consuming issue of adap- 
tive grid-refinement. However, the method has its weaknesses 
compared to grid-based methods. For example, its conver- 
gence behavior is mathematically difficult to assess and the 
method has problems reproducing certain types of dynami- 
cal instabilities ( [Agertz et al. \ |2007| l. The algorithmic sim- 
plicity of the method and the high flexibility due to its La- 
grangian nature, however, usually outweigh these drawbacks 
and SPH remains one of the numerical workhorses of current 
star-formation studies. There are various implementation of 
the method. Examples of very popular codes are GADGET 



(|Sprin gel 2005 ; Springel et al. 20 l|l, GASOLINE ( fWadsley 



et al.\ |2004) , MAGMA (Rosswog and Pricej |2007| l, and the 



various decendents of the SPH program originally developed 
by Be nz ( |Bate et al.\ [1995) [Bate and BuAertl [T9971 [Benz} 



3. Sink Particles as Subgrid-Scale IVIodels for Protostars 

A fundamental problem for modeling protostellar collapse 
and star formation are the enormous density contrasts that 
need to be covered. Regions of high density require small 
cell sizes in grid-based methods (Truelove et al. 1997[), or 



equivalently, small-particle masses in SPH calculations ( Bate 
[and BurkertI |1997t [Whitworthj |1998| l. In order to guarantee 
stability, every numerical scheme must resolve the traversal 
of sound waves across the minimum resolution element, i.e. 
either across one cell or across the smoothing kernel of indi- 
vidual SPH particles. This is the so called Courant Friedrich 
Lewy criterion. It causes the time integration stepsize to get 
smaller and smaller as the density increases. As a conse- 
quence, a computation virtually grinds to a halt during gravi- 
tational collapse. 

When modeling the build-up of entire clusters of stars, or 
even following the accretion of the bulk of a core's mass onto 
a single star, this problem clearly needs to be overcome. One 
way out is to introduce sink particles. Once the very center of 
a collapsing cloud cores exceeds a certain density threshold 
(usually several thousand times the mean density, or using a 
threshold based on the Jeans mass) it is replaced by one single 
particle which inherits the combined masses, linear and an- 
gular momentum of the volume it replaces and which has the 
ability to accrete further gas from the infalling envelope. This 
permits to follow the dynamical evolution of the system over 
many global free-fall times, however, at the cost of not be- 
ing able to resolve the evolution at densities above the thresh- 
old value. In a sense, sink particles introduce "inner bound- 
aries" to the computational domain. They have been success- 
fully introduced to grid-based (Krumholz et flZ.','2004 1 as well 



as particle-based methods (Bate et al. 1995; Jappsen et al. 
12005] !. ~ 

Each sink particle defines a control volume with a fixed 
radius. It lies typically between a few and a few hundred 
astronomical units, AU, depending on the specific goals of 
the calculation. For comparison, the radius of Earth's or- 
bit per definition is exactly 1 AU. In most cases the sink ra- 
dius is chosen such that the Jeans scale below the thresh- 
old density is sufficiently resolved ( |Bate and Burkertj |1997[ 
[Truelove et a/.[[T997] l. There are, however, implementations 
where sink particles have radii equivalent to one cell only. Be- 
cause the interior of the control volume is not accessible, the 
physical interpretation is often very difficult and subject to 
debate. Usually sink particles are thought to represent indi- 
vidual protostars or dense binary systems. This is supported 
by detailed one-dimensional implicit radiation hydrodynamic 
calculations which demonstrate that a protostar will build up 
in the very center of the control volume about 10^ yr after 
sink creation (Wucht erl and Klessen| |2001| l which will swal- 
low most of the infalling material. 

Protostellar collapse is accompanied by a substantial loss of 
specific angular momentum, even in the absence of magnetic 
braking (,Fisher,,2004, Jappsen and Klessen , 2004) . Still, most 
of the matter that falls in will assemble in a protostellar disk. It 
is then transported inward by torques from magnetorotational 



and possibly gravitational instabihties (Bodenheimer 1995 
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Gammie' '2001', 'Kratter and Matzner', '2006', 'Kratter et al 
2008; Laughlin and Bodenheimer, 1994; Lin and Papaloizou 



1996 



et al. 



jLodato and Rice 2005; Papaloizou and Lin| |1995) Shu| 
)l990^ ). With typical disk sizes of order of several hun- 



dred AU in simulations of the formation of star clusters, the 
control volume fully encloses both star and disk. Even in 
higher resolution calculations that focus on single cores, the 
control volume contains the inner part of the accretion disk. If 
low angular momentum material is accreted, the disk is stable 
and most of the material ends up in the central star. In this 
case, the disk simply acts as a buffer and smooths eventual 
accretion spikes. It will not delay or prevent the mass growth 
of the central star by much. However, if material that falls 
into the control volume carries large specific angular momen- 
tum, then the mass load onto the disk is likely to be faster 
than inward transport. The disk grows large and may become 
gravitationally unstable and fragment. This may lead to the 

formation of a binary or higher- order multiple ( Bodenheimer| [and Bonnelll. [20081 |McKee and Ostriker||2667HZinnecker and 



these stars radiate prodigiously via Kelvin-Helmholtz contrac- 
tion and then nuclear burning. The first effect of the radiation 
they produce is on their immediate environs. Their starlight 
is absorbed by dust grains suspended in the circumstellar gas, 
exerting a pressure that opposes gravity. Second, as the ra- 
diation diffuses out of the dusty gas clouds around a massive 
star it heats the gas. This affects the process of fragmentation, 
and thus plays a role in determining the stellar mass function 
for all stars born in strongly irradiated environments. Third, 
once massive stars contract onto the main sequence, they be- 
come significant sources of ultraviolet radiation, which can 
dissociate molecules, ionize atoms, and drive strong shocks 
throughout the star-forming cloud. These processes can both 
inhibit and promote star formation. 

Radiation Pressure in Massive Star Formation. The first 
effect is perhaps the best studied, and has been the subject 
of several recent reviews ( Beuther et al.\ |2007[ Krumholz 



et al. 2000 Kratter et al. 2008 1. Indeed a initial binary frac 



tion of almost 100% is consistent with observations of star 
clusters ( Kroupa| |1995| l. To some degree this can be taken 
into account by introducing an appropriate scaling factor. In 
a cluster environment the protostellar disk may be truncated 



by tidal interactions and loose matter (Adams et al. 2006 



Clarke et al. 20001. The importance of this effect depends 



strongly on the stellar density of the cluster and its dynam- 
ical evolution. Further uncertainty stems from the possible 
formation of O or B stars in the stellar cluster. Their intense 
UV radiation will trigger evaporation and gas removal, again 
limiting the fraction of sink particle mass that turns into stars. 
Similar holds for stellar winds and outflows. These feedback 



processes are discussed in Section IV below. 



IV. THE IMPORTANCE OF FEEDBACK 

A. Feedback Processes 

Most star formation simulations to date have neglected the 
effects of feedback on the star formation process. While this 
is computationally simpler, it is clearly not physically correct, 
and its omission leads to a number of obvious differences be- 
tween simulations and observations. For example, simulations 
without feedback produce star formation that is too rapid and 
efficient ( Krumholz et al. 2007a| l and significantly overpro- 
duce brown dwarfs compared to observations ( [Offner et al.\ 
[2008'). To make progress the next generation of simulations 
will have to remedy this omission. 



1 . Radiation Feedbacl< 

We begin our consideration of feedback processes by exam- 
ining the effects of radiation from young stars. It is convenient 
to distinguish three distinct types of radiative feedback on star 
formation. The dominant sources of radiation in forming star 
clusters are the young massive stars, which begin their lives 
accreting rapidly, producing high accretion luminosities from 
the initial infall onto their surfaces. Later on in their formation 



Yorke 2007j , so we skip over it relatively quickly. As early 



as the 1970s researchers considering the formation of mas 
sive stars realized a fundamental problem. The largest stars in 
nearby galaxies have masses ~ 100—150 Mq (Bonanos etal. 
2004} |Figer| [2005] |Rauw et al.\ |2005| l, and for these stars ra 



diation pressure on electrons within the star is the dominant 
support mechanism. In effect, these stars are at their inter- 
nal Eddington limit. However, the Thompson cross-section is 
smaller than the cross-section of dusty gas to stellar radiation 
reprocessed into the infrared by an order of magnitude. Thus 
if a massive star is at its Eddington limit with respect the ion- 
ized, dust-free gas in its interior, it must exceed the Eddington 
limit by an order of magnitude with respect to the dusty gas 
found in molecular clouds. How then is it possible for dusty 
gas to accrete and form a massive star, since the outward radi- 
ation force on the accreting material should be significantly 



stronger than th e pull of gravity (|Kahn| 1974[ Larson and 



Starriield|[T97Tl|Wolfire and Cassinelli||1987[|Zinnecker and 



'^rkel|2007^ 

Analytic treatments of the problem suggest that the solution 
lies in the non-sphericity of the accretion process: if the dusty 
gas around a protostar is sufficiently opaque, it can coUimate 
the radiation, reducing the radiation force over some fraction 
of the solid angle to the point where gravity is stronger and ac- 
cretion can occur ( Jijina and Adams 1996; Krumholz et al. 



2005b[ |Nakano| |1989,. Nakano et al, 1995 ). Simulations ap- 
pear to bear out this solution, at least preliminarily. Hydro- 
dynamic simulations in two dimensions using a flux-limited 
diffusion approach (see below) are able to form stars up to 



40 Mq before radiation pressure reverses infall (Yorke and 
|Sonnhalter| |2002| ), while three-dimensional simulations show 
no signs of a limit on the upper masses of stars imposed by 



radiation pressure ( Krumholz et al. \ |2005a[ |2009). In both 
the 2-D and 3-D cases, radiation is strongly beamed toward 
the poles of an accretion disk, allowing gas to accrete through 
parts of the equatorial plane shielded by the disk. In 3-D, 
this self-shielding is further enhanced by the organization of 
the gas into opaque, dense filaments, while radiation escapes 
through optically thin channels. This effect appears to allow 
the formation of stars with no clear upper mass limit. 
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Radiation Heating and the IMF. The second radiative ef- 
fect is heating of the gas, with the resulting modification of 
the initial mass function. Increasing the gas temperature sup- 
presses fragmentation, and the observed overproduction of 
brown dwarfs in isothermal simulations ( |Bate et aT] [2003 | l 
is at least in part due to their omission of radiative feedback 
( [Matzner and Levin| |2005| l. Similarly, radiative feedback is 
a strong candidate solution to another mystery about mas- 
sive star formation: why would ^ 100 Mq of gas, a mass 
that represents tens to hundreds of Jeans masses at the typi- 
cal densities and temperatures of a molecular cloud that has 
not yet begun to collapse, ever collapse coherently rather than 
fragmenting into many objects (Bonnell et al. 2006[ Dobbs 
\et al.\ |2005| l? A possible answer is that the accretion lumi- 
nosity produced by the collapse of a dense core in a massive 
star-forming region is sufficient to suppress a high level of 
fragmentation, converting a collapse that might have produces 
~ 100 small stars into one that produces only a few massive 
ones ( [Krumholzl [20061 [Rmmholz and McKee[[2008] l. How- 
ever, the overall importance of this process and its details are 
subject to ongoing debate, and clearly more work is required 
on this important subject. 

To date there is only one published simulation of the ef- 



fect of radiative heating on fragmentation (Krumholz et al. 



2007a I, and it confirms the analytically-predicted outcome. 



Radiative heating reduces the amount of fragmentation that 
occurs during the collapse of massive pre-stellar cores. Figure 
|9] shows an example of this effect, comparing two simulations 
that are identical in every respect except that one is done with 
radiative transfer to one without it. However these simula- 
tions only studied single, isolated cores, and thus tell us little 
about the effect of radiative feedback on the IMF in star clus- 
ters. It remains unclear how radiative feedback in star clusters 
shapes the IMF. However it seems clear that, in the dense, 
optic ally- thick environments where clusters form, any results 
derived from simulations using the isothermal or optically- 
thin cooling approximations will increasingly diverge from 
reality once high-mass stars begin to form. 

High Energy Radiation and Star Formation Efficiency. The 
third form of radiative feedback from massive stars is high 
energy radiation that is capable of dissociating hydrogen 
molecules (photon energies above 11 eV) and ionizing hy- 
drogen atoms (photon energies above 13.6 eV). The former 
creates a photodissociation region (PDR), a volume of mixed 
atomic and molecular gas at temperatures of hundreds of 
Kelvin, too warm to form stars. The latter rapidly heats the 
gas around a massive star to ^ 10"* K and raises its sound 
speed to ~ 10 km s^^, forming a structure known as an H II 
region. Except in the case of stars with very weak ionizing 
fluxes, or in environments where the magnetic field strongly 
confines the ionized region, the shock front generated by an 
expanding H II region generally overruns the PDR created by 
dissociating radiation and traps it between the ionization front 
and the shock front. For the purpose of molecular cloud dy- 
namics, therefore, ionizing radiation is usually the more im- 
portant effect ( [Rrumholz et a/.|[2007c| l. 

Unlike radiative heating and radiation pressure, dissociat- 
ing and ionizing feedback do not become significant until 



fairly late in the star formation process. Early on rapid ac 
cretion swells massive stars to radii of ~ 100 Rq ( Hosokawa 



|and Omukai||2008")|Yorke and Bodenheimer|[2008[ l, and these 
large radii lead to low surface temperatures, reducing the frac- 
tion of a massive star's power that emerges at energies above 
13.6 eV. Moreover, even the full main sequence ionizing lu- 
minosity from a massive star will not escape from the stel- 
lar vicinity if accretion onto the star is sufficiently rapid and 



covers enough of the stellar surface (Keto 2002 2003 Tan 



.and McKee^ 2003, Walmsley , 1995) ). In this quasi-spherical 
case, the H II region is kept from expanding and the Strom- 
gren radius is small. However, the situation may change once 
protostellar outflows are taken into account (see next Section 
|IV.A.2[ ). These effectively remove high-density material along 
the rotational axis of the system. This may lead to an H II re- 
gion that escapes along the outflow axis while remaining con- 
fined in the equatorial direction ( |Tan and McKee| |2003| , or 
it may allow the H II region the break free entirely from its 
parent core (Keto 2007| l. 

Once ionization does begin to break out of a massive pro- 
tostellar core, it is likely to be the most significant of the 
three types of radiative feedback. Since 10 km s^^ is much 
greater than the escape speed from a molecular cloud under 
Milky Way conditions, ionized gas escapes from star-forming 
clouds into the ISM, reducing the amount of mass available 
for star formation and unbinding molecular clouds (McKeej 
land WilUamsl [19971 fWilliams and McKee, .1997j . Further- 
more, since 10 km s~^ is much larger than the sound speed 
in the non-ionized molecular gas, once they form Hll regions 
expand dynamically, driving shocks into the neutral material. 
Analytic models suggest that this can both promote star for- 
mation, by sweeping up gas into sheets that subsequently frag- 



ment by gravitational instability (Eknegreen and Lada 1977 



Whitworth et aL \ [1994||, a nd inhibit it, by driving turbulent 



motions ( [Krumholz etaL\ [2006[ [Matzner| [2002| l. Clearly 
more work is required to determine which effect dominates. 

Simulations of these processes are still quite primitive, and 
most have focused on small molecular clouds that are already 
in a process of free-fall collapse when the simulation begins. 
Within this limited context simulations have produced a num- 
ber of qualitative conclusions. First, single ionizing sources at 
the molecular cloud centers do not easily unbind those clouds, 
even if they deposit an amount of energy larger than the cloud 



binding energy (Dale et a/.||2005 ; Mac Low et al. 20071. This 
is because in a cloud with a pre-existing density structures, 
most of the energy is deposited in low-density gas that freely 
escapes from the cloud, while the higher density material is 
largely unaffected. Thus, in this context the effects of ion- 
ization on reducing the star formation efficiency are modest. 
Second, ionization can drive significant velocity dispersions 



in neutral gas, possibly generating turbulence ( jMellema et al. 



[2006[ ). Third, ionization does sweep up material and promote 
coflapse, but numerical simulations indicate that this effect 
may again be modest (Dale et al. 2007a'bi'Mac Low et al. 



,2007 ) . Much of the swept up gas in these calculations was al- 
ready on its way to star formation due to gravity alone, and the 
compression produced by an H II region shock only modifies 
this slightly. 
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FIG. 9 A comparison of two simulations witli identical initial conditions and evolution times, one including radiative transfer {left panel) and 
one done without it (right panel). Stars are indicated by plus signs. The simulation without radiative transfer forms a factor of ~ 4 more stars 
than the one including it, and has significantly less mass in its gaseous disk. (Images adapted from \Krumholz et aL] ( |2007fl| >J. 



This work is only a beginning, and many questions remain. 
First, none of the simulations to date have included multi- 
ple ionizing sources that are simultaneously active, so that 
interactions between expanding H II region shells can pro- 
mote both star formation and turbulence. Since massive stars 
form in clusters, however, multiple sources should be the rule 
rather than the exception. Second, the simulations have for 
the most part focused on small, tightly-bound proto-cluster 
gas clouds being ionized by rather small ionizing luminosities 
corresponding to single stars, rather than larger, lower den- 
sity, more loosely bound molecular clouds subjected to the 
ionizing flux of an entire star cluster. The effects of ioniz- 
ing radiation may be greater in the latter case than in the for- 
mer. Finally, only one simulation of ionizing radiation feed- 



back to date has included magnetic fields CKrumholz et al. 



[2007c ), and only then in a very idealized context. Since mag 
netic fields can tie together high- and low-density regions of a 
cloud, they may significantly increase the effects of ionization 
feedback. 



2. Protostellar Outflow Feedback 

Outflows from young stars provide another significant 
source of feedback on local scales in star-forming regions. 
During the process of accretion onto young stars about ^ 10% 
of the gas that reaches the inner accretion disk is ejected into 
a collimated wind that is launched at a speed comparable to 
the Keplerian speed close to the stellar surface. Theoretical 
predictions ([Anderson et a l.\ "SOUS'; 'Konigl and Pudritz', '2000', 
[PelletierandPudritz|[1992„Pudritz„2003„Pudi-itz efa/.„2007, 



'Shu et al' '2000 1 and observational data ( [Bontemps et al.\ 
,1996; Richer et al. 2000^ agree very well on this value. An- 
other quantity that is well constrained by observations is the 
net momentum flux of the material entrained in the outflow. 
It is typically p ~ 0.3Mvk, where M is the accretion rate 
onto the star plus disk and vk is the Keplerian velocity at the 
stellar surface (Bontemps et al. 1996 Matzner 2002[ Richer 
\et al. \ [2000p . Outflow momentum flux correlates well with 
source luminosity across a very wide range in luminosity L, 
suggesting that all protostars show a common wind launching 
mechanism independent of mass ( Wu et aL\\2Q04^ . Since the 
wind momentum flux is much greater than L/c, this mecha- 
nism is almost certainly hydromagnetic rather than radiative in 



nature. The two primary theories for this are the x-wind ( Shu 
etal."2000) and the disk wind ([Konigl and Pudritz|[2000||Pir 



dritz 2003 ; Pudritz et al. \ [2007| l. Common to both models 



is the idea that matter gets loaded onto magnetic field lines 
and then accelerated outwards by centrifugal forces. From the 
standpoint of feedback on scales large compared to the ac- 
cretion disk around the source, the details of how the wind 
is launched matter little. All magneto-centrifugal winds ap- 
proach the same distribution of momentum flux per unit angle 



at large distances from the launching region (Matzner an d Mc 
|Kee[[1999] ). On larger scales, on which the outflow interacts 
with ambient material in the core, the opening angle varies 
depending on mass and age of the protostar. Outflows from 
low-mass stars appear quite well collimated, and remain so up 
to roughly B stars ([Arce et aL\ [2007 [ [Beuther and Shepherd[ 



2005 Richer et al. 2000 1. The opening angles of O star out- 
flows are wider, but it is unclear if this widening is an inherent 
property of the outflow or a result of the interaction between 
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the outflow and the ambient gas. 

One important difference between outflow and radiation 
feedback is that outflow feedback is more democratic. The 
most massive stars in a cluster dominate its radiative output, 
because (except at the very highest stellar masses) luminos- 
ity is a very strong function of mass, and ionizing luminosity 



an even stronger one (Kippenhahn and Weigert 1994) ). The 
dependence of luminosity on mass is strong enough to over- 
come the relative dearth of massive stars compared to low- 
mass ones. For outflows the reverse is true. The total mass 
accretion rate onto all the stars in a cluster is necessarily dom- 
inated by the low-mass stars, since they comprise the bulk of 
the stellar mass once star formation is complete. The Keple- 
rian velocity at a star's surface varies as ^JM/R, where M 
and R are the star's mass and radius, and this ratio is only a 
very weak function of mass for main sequence stars. Thus, 
we might expect low-mass stars near the peak of the IMF to 
dominate outflow feedback. This simple analysis neglects the 
effect that more massive stars have shorter Kelvin-Helmholtz 
times and thus reach smaller radii more rapidly than low-mass 
stars, giving them larger Keplerian velocities at earlier times. 
Including this effect in a more careful analysis suggests that 
each logarithmic bin in mass contributes roughly equally to 
the total amount of momentum injected into a cluster ( ,Tan and, 
|McKee|[2002| l. This has two important consequences: first, it 
means that outflow feedback can be important even in small 
clumps that do not form massive stars. Second, it means that 
simulations of outflow feedback cannot focus exclusively on 
the most massive stars, but must instead consider all stars as 
sources. 

Outflows can influence their immediate surroundings as 
well as the cluster in which they form. On small scales, they 
reduce the star formation efficiency by removing mass from a 
collapsing core, both directly and via material that the outflow 
entrains as it escapes the core. Analytic estimates su ggest that 
this process removes 25 — 75% of the mass in a core (Matzner 



[and McKee| |2000| , but this is highly uncertain since no sim- 
ulations of the collapse of individual cores with outflows that 
are capable of evaluating this estimate have been published. 

Outflows could also modify the star formation process by 
removing mass from protocluster gas clumps and possibly by 



driving turbulence within them ( 


Matzner 


2007 


Matzner and 


IMcKeel |2000l |McKee| |1989| |Norman anc 


Silk 


11980 1. How- 



ever, it must be noted that this process cannot be the main 
driver of turbulence on global molecular cloud scales as out- 
flows typically have short length scales. Instead this turbu- 
lence is probably driven on scales comparable to or larger 
than typical cloud sizes ( [Mac Low an d Ossenkopf 2000j 
[Ossenkopf efal] |2008a|b| [OssenkopTand Mac Low, |2002l l. 
One possible candidate for the origin of this large-scale tur- 
bulence is convergent flows in the galactic disk (see Section 



II.B.2I either driven by gravitational instability in the disk 
( |Dobbs et al.\ [2008] |Li et a/.| [2005| [2006|, by collisions be- 
tween molecular clouds ( Tasker and Tanj 2009 1, or caused by 



supernova explosions (de Avillez and Breitschwerdt 2007 



[Mac Low et al.\ |2005| l. Another candidate is giant H II re 
gions created by clustered star formation within clouds, which 
have size scales comparable to entire GMCs ( jKrumholz et al. 




FIG. 10 Evolution of the total kinetic energy (upper line) and grav- 
itational energy (lower line) as a function of time in a simulation of 
the formation of a star cluster including protostellar outflow feed- 
back. Energies are normalized to the initial kinetic energy in the 
simulation, and times to the gravitational (or free-fall) time. (Plot 
taken from \NaIcamura and Li\ ^2007\ ) 



|2006t [Matzner! [2002| l. In addition, there seems to be no dif- 
ference between the measured turbulence content of cloud 
clumps that are still in the so-called dark phase, i.e. before 
star formation has set in, and cloud clumps that are already 



actively building up stars in their interior (Ossenkopf et al. 
|2001| ). This indicates that, at least at birth, star-forming re- 
gions must have turbulent motions that were imprinted as 
part of the formation process. Conversely, however, obser- 
vations show that high column density star-forming clumps 
within GMCs lie above the linewidth-size relation observed 



for GMCs as a whole (Heyeref al. 2008 Plume et al. 1997 



Shirley et al. 2003 | l. This suggests that their turbulence can- 



not be supplied from large scales motions within the parent 
GMCs. Either these regions are powered by gravitational col- 
lapse, in a scaled-down version of the scenario described in 



Section 11. B. 2 or they are driven by internal sources. Out- 
flows are a natural candidate for this, and the deviation from 
a simple powerlaw linewidth-size relation predicted by ana- 
lytic models appears to be consistent with what is observed 
(Matznerl [20071 ). 

Simulations of protostellar outflows to date fall into two 
categories. Local simulations focus on the interaction of a sin- 
gle outflow with an ambient medium at high resolution, while 
larger-scale simulations follow an entire gas clump and star 
cluster including multiple outflows, but at significantly lower 
resolution. Local simulations attempt to understand the driv- 
ing of turbulence by single outflows in detail. However, in- 
terpretation of these results is difficult, since there is no sim- 
ple way to separate "turbulence" from the coherent motion 
caused by a single outflow. Different authors analyzing sim- 
ulations in different ways have come to opposite conclusions, 
with some arguing that outflows cannot drive supersonic tur- 



can ( [Cunningham ef g/. [2008 



bulence (Banerjee et al.\ 2007 i, while others conclude that it 



2006| l. 



Global simulations of outflow feedback generally find that 
it has strong effects on the cluster formation process. In 
the absence of energy input, simulations of cluster-forming 
gas clumps find that any turbulence initially present decays 
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rapidly, leading to a global collapse in which an appreciable 
fraction of the mass is converted into stars within a few dy- 
namical timesjBate et al. 2003 ; Bonnell et al. 



2003] [KL 



essen 



land Burkert] [2000, |2001t |Klessen et al.\ |199 8i Simulations 
that include outflow feedback found that outflows can change 
this picture. They eject mass from the densest and most ac- 
tively star-forming parts of a cluster, reducing the star forma- 
tion rate, while at the same time injecting enough energy to 
slow down overall collapse and maintain a constant level of 
turbulent motions {Li and Nakamura] [2006 [ [Nakamura and] 
[LH |2007 |2008| I. As a result, the star formation rate drops 
to < 10% of the mass being converted into stars per free- 
fall time, and rather than undergoing a runaway collapse the 
clump reaches a slowly evolving quasi-equilibrium state. Fig- 
ure [To] illustrates this effect in a simulation of the formation 
of a star cluster including protostellar outflow feedback. At 
the start of the calculation the kinetic energy falls as the initial 
turbulent velocity field decays. Consequently the cloud con- 
tracts and at about half of a free-fall time stars start to form 
and drive outflows. This energy input changes the subsequent 
evolution, and the cloud's global contraction is halted or at 
least significantly retarded. A definite answer would require 
to follow the dynamics over a longer period in time. 

All the simulations published to date have significant lim- 
its. With only one exception they treat the wind as an in- 
stantaneous explosion, rather than a continuous beam injected 
over ~ 10^ yr as we observe. The individual explosion spikes 
are clearly visible in Figure 10 The current studies are also 



characterized by low numerical resolution (128^ cells), which 
makes the energy injected more space filling, which is one 
of the main characteristics of interstellar turbulence. In ad- 
dition, the calculations are performed in a periodic box, so 
energy cannot leave the star-forming cloud. In reality some 
very strong pencil-beam outflows escape their parent clumps 
and cover distances of a few parsec ( B ally^ [2007| [Stanke etal.\ 
2002| l. The one simulation published thus far that does include 



time history of accretion and better resolution only considers 
the effects of outflows from stars larger than 10 Af© ([Dale and 



[Bonnell|[2008[ ), thereby neglecting the majority of the outflow 
power. Clearly the problem of outflow feedback and how it af- 
fects star formation is in need of further study. 



3. Other Types of Feedback 

Although radiation and protostellar outflows are thought to 
be the dominant feedback processes in star formation, two 
other mechanism are worthy of brief discussion: supernovae, 
and winds ejected by stars on the main sequence and post- 
main sequence. In terms of sheer energetics, it might seem 
odd to ignore supernovae as a major source of feedback. How- 
ever, two compensating effects reduce their role in regulating 
star forming clouds. The first is timescales. Even the most 
massive stars do not explode as supernovae until 3 — 4 Myr 



after formation ( Parravano et al. 2003j), and this is compara- 
ble to or longer than the formation time of star clusters. Thus 
supernovae come too late to affect the formation of individual 
star clusters, although they may be able to affect theii^ parent 



giant molecular clouds, which have longer lifetimes. 

A second effect, however, mitigates the impact of super- 
novae on GMC scales as well. Supernovae occur only after 
Hll regions and stellar winds have carved large cavities of hot, 
ionized gas around the massive stars that produce them. If a 
supernova occurs while this bubble is still embedded within 
its parent cloud, much of its energy is radiated away while 
the blast wave is still confined to the bubble. Simulations find 
that, as a result, the mass that is removed from the cloud by 
a supernova plus ionization is typically only ~ 10% larger 



than that removed by ionization alone ( 


Krumholz et al. 


2006[ 


Matzner 2002 Tenorio-Tagle ef a/. 1985| 


Yorke et al. 


19891. 



If, on the other hand, the bubble of hot gas created by ion- 
ization has broken out of a massive star's parentel cloud by 
the time the star explodes, both distance and an impedance 
mismatch make it difficult to deliver much of the supernova 
energy to the cloud. In a few Myr the expanding Hll region 
around a massive star cluster can push back the parent GMC 
by ^ 10 pc or more from the site of the supernova, which then 
occurs in an ionized medium whose density is 2 — 4 orders of 
magnitude lower than that of the cloud. Both the distance 
and the large density jump serve to shield a molecular cloud 
from the effects of a supernova, so that very little of the super- 
nova energy is deposited in the molecular cloud. This effect 
means that, even on GMC scales, supernovae are unlikely to 
be the dominant feedback mechanism. It is important to note, 
however, that the energy that is not deposited in the molec- 
ular cloud itself does nevertheless affect the remainder of the 
ISM on galactic scales. Consequently supernovae are likely to 
dominate the energetics of the ISM on large scales ( |Mac Low] 
[and Klessen||2004| l. 

Main sequence winds are also thought to be subdomi- 
nant as feedback mechanisms due to the effects of ionization 
(|Matzner| [2002[ [McKee et oT] [T984). Stellar winds initially 



expand into a bubble of ionized gas created by the ioniza- 
tion from a massive star, and they create a radiative shock 
within that bubble where much of the wind energy is dissi- 
pated. Only after the stellar wind shock catches up to the 
ionization-created shock can stellar winds begin to provide 
feedback to the parent cloud. Even then the increase in to- 
tal kinetic energy in the shock is modest for stars up to at least 
35 Mq ( [Freyer et al. 2006 1, and even for 60 Mq stars is only 
of order unity ( Freyer et al. 2003 i. 



B. Modeling Feedback 

1 . Numerical Methods for Non-Ionizing Radiation Feedback 

Stars emit the bulk of their radiation in the visible part of the 
spectrum, but the dusty clouds in which stars form are gener- 
ally very opaque to visible light until late in the star formation 
process, when most of the gas has already been accreted or 
dispersed. As a result, direct stellar radiation tends to be ab- 
sorbed by dust and reprocessed into the infrared close to the 
star that emits it, and modeling the resulting diffuse infrared 
radiation field is the primary goal of most numerical methods 
for simulating non-ionizing radiation feedback. 
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Focusing on the diffuse infrared radiation field simplifies 
the radiative transfer problem considerably, since the primary 
opacity source at infrared wavelengths is dust rather than 
atomic or molecular lines, and because in the IR scattering 
is negligible compared to absorption ([Rybicki and Lightman 



|T979|. Even with these simplifications, though, it is possible 
to solve the full equations of (magneto-)hydrodynamics plus 
the equation of radiative transfer for this problem only in one 
dimension (Larson 1969 Masunaga and Inutsuka' 2000; Ma-| 
|sunaga eFai, 1998j . Such an approach is unfortunately too 
computationally expensive to be feasible in three or even two 
dimensions. Instead, one must simplify the problem even fur- 
ther 

One approach is simply to modify the standard optically- 
thin cooling curve used in simulations without feedback by 
using an approximation to estimate the optical depth and re- 
duce the cooling rate appropriately ( Banerjee etal. 2006[ Sta- 
[matellos et al.\ pOOT aV In this case one need not to solve a 
radiative transfer problem at all. This is an advantage, since it 
means that the radiation step is has nearly zero computational 
cost, but it is also a limitation. Because it lacks a treatment 
of radiative transfer, this approach allows gas to heat up due 
to adiabatic compression, but not because it is being illumi- 
nated by an external radiation source. In particular, in this 
approach there is no way for stars to heat gas. Since stellar ra- 
diation provides significantly more energy than gravitational 
compression once the first collapsed objects form (Krumholz 



20061 IKrumholz et al. 2007aj [Krumholz and McKee, |2008 



Yorke and Sonnhaltert|2002t|Zinnecker and Yorke(|2007l l, this 



technique is only suitable for simulating star formation up to 
the point when the first parcels of gas collapse to stars. 

The most common and simplest approach that can go past 
first collapse and follow either accretion or subsequent star 
formation, and the only one used so far in "production" simu- 
lations ("Krumholz et flZ.','2005a 



2006 , Yorke and Bodenheimer 



2007a 



1999^ 



JWhitehouse and Bate 
Yorke and Sonnhalter 



2002|l, i s the flux-limited diffusion approximation ( |Alme anc 



Wilsonj |1973[ |Levermore and Pomraning"] |1981| l. The un- 



derlying physical idea is simple: in an optically thick envi- 
ronment like the dusty clouds in which stars form, radiation 
diffuses through the gas like heat, and the radiative flux F 
obeys Pick's Law: F = ~cV E / {"inp), where E is the radia- 
tion energy density, p is the gas density, and k, is the specific 
opacity, with units of area divided by mass. This is the stan- 
dard diffusion approximation, and it can be made either in a 
gray form by integrating E and F over all frequencies, or in a 
multi-group form in which one divides the spectrum into some 
number of intervals in frequency and computes a separate en- 
ergy density and flux for each interval ( [Shestakov and Offner[ 
[2^ [Yorke and Sonnhalter] [2002l i. 

The pure diffusion approximation encounters a problem 
when the opacity is low, since if Kp is sufficiently small the 
flux can exceed cE, violating the constraints of special rel- 
ativity. The flux-limiting approach is to solve this problem 
by modifying the law for the flux to F = —Xc\/E/{kp), 
where A is the flux limiter, a dimensionless function of E and 
Kp that has the properties A — > 1/3 when the gas is opti- 
cally thick, and A KpE/\VE\ when it is optically thin. 



This limiting behavior ensures that flux approaches the cor- 
rect Pick's Law value when the optical depth is high, and cor- 
rectly reaches a maximum magnitude of cE at low optical 
depth. Many functional forms are possible for A. The most 
commonly-used one is the Levermore & Pomraning limiter 
A = R-^{cothR - R-^), widi R = \yE\/ {KpE) ( |Lever- 
|more| [T984[ [Levermore and Pomraning||198l] l. 

Given a formula for computing the radiation flux in terms of 
the radiation energy density, it is possible to drop all moments 
of the equation of radiative transfer except the zeroth one, so 
that the set of equations to be solved consists of the standard 
equations of HD or MHD, with some added terms describing 
the interaction of radiation with the gas, plus one additional 
equation for the radiation energy density. One treats feedback 
from stars in this formulation simply by adding it as a source 
term or a boundary condition in the radiation energy equation. 
The resulting set of equations may be written using either a 
comoving ("Boss and Myhill' 799?; 'Hayes et al.\ '2006', "Stoiie 
iri992 



\etaL\ \1992: Tscharnuter, 1987; Whitehouse and Bate 2004J 
Whitehouse et al. 20 0j]l or a m ixed-frame (Howell and Gre^ 
nough 2003 , ^Krumholz et al. 2007b Shestakov and Offne^ 
2008 ) formulation. The former approach is more suited to im- 
plementation in a code that is either Lagrangean, such as SPH, 
or based on van Leer advection ( |van Leer] |1977j , but has the 
disadvantage that the equations are not expUcitly conservative, 
and so the resulting codes cannot precisely conserve energy. 
The mixed-frame equations, on the other hand, are explicitly 
conservative, which makes them preferable for codes based 
on a conservative update, particularly those involving adap- 
tive mesh refinement. In either form the equations are still 
significantly more expensive to solve than the corresponding 
non-radiative ones, since codes must handle radiative diffu- 
sion implicitly in order to avoid severe constraints on the time 
step, but solutions are within the reach of modern supercom- 
puter simulations. 

Although it is the tool of choice for star formation simula- 
tions at present, the pure flux-limited diffusion approximation 
does have some important limitations. Diffusion methods do 
not correctly represent shadowing effects which appear in sys- 
tems that are optically thin or nearly so, nor can they model 
direct stellar radiation before it is absorbed and re-radiated 
isotropically. Pure diffusion also assumes that the gas and 
dust are thermally well-coupled. While this approximation is 
a good one at densities 10^ cm"'^ or more, it may fail at 
lower densities. Diffusion also neglects cooling via molecular 
line emission, which can also be important at lower densities 
( |Genzel|[T99T] l. 

The literature contains a variety of numerical techniques 
to address these shortfalls. One can handle imperfect dust- 
gas coupling by explicitly including it in the iterative radia- 
tive transfer update ( [Yorke and Bodenheimer 1999 1. To han- 
dle direct stellar radiation or molecular cooling as well as the 
diffuse IR field, one can use a hybrid approach that com- 



bines a diffusion step with a ray-tracing step (Murray et al. 



19941 or an optically-thin cooling step. To correctly model 



shadowing, one can use a more sophisticated radiative trans- 
fer method than diffusion, such as Monte Carlo, ray-tracing 
([Heinemann et al.\ [2006[l, variable tensor Eddington factor 
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(VTEF) (Hayes and Norman 2003 | l, or 5„ transport (Livne 



et al. 20041. However, with the exception of the dust-gas et al. 2007; Mellema et al. 



coupUng method, none of these techniques have thus far been 
used in any "production" simulations of star formation. In 
some cases this is simply a matter of the necessary techniques 
not yet having been implemented into the codes most com- 
monly used for star formation studies. These techniques, such 
as two-step approaches for the diffuse IR field and direct stel- 
lar fields and line radiation, are likely to appear in production 
simulations in the next few years. In other cases, however, the 
limitation is one of computational expense. For example the 
VTEF and Sn methods have thus far only been used in two- 
dimensional calculations, simply because in three dimensions 
they have thus far proven prohibitively expensive. Remedying 
these problems will require significant advances in radiative 
transfer methodology to solve. 



2. Numerical Methods for Ionizing Radiation 

In comparison to non-ionizing radiation, handling ionizing 
radiation is conceptually more straightforward. Rather than 
a diffuse field arising from the repeated reprocessing of stel- 
lar radiation by dust grains, the ionizing radiation field in a 
star-forming region consists mostly of photons directly emit- 
ted from a stellar surface. Only in the outer parts of low- 
density ionized regions with sharp density gradients does re- 
processed radiation make up a significant part of the photon 
field ( [Ritzerveld] |2005[ ). The dominance of a relatively small 
number of point sources of radiation translates into a much 
simpler computational problem. By far the most common ap- 
proach for solving it is to adopt the on-the-spot approxima- 
tion ( jOsterbrock, 19^9 ), in which one assumes that recombi- 
nations of ionized atoms into the ground state yield ionizing 
photons that are re-absorbed immediately near the point of 
emission. One therefore ignores photons emitted by recom- 
bining atoms entirely, and one solves the transfer equation 
along rays from the emitting stellar source, balancing recom- 
binations into excited states against ionizations along each ray. 

Within this over-arching framework, there are a variety of 
subtleties about how one draws the rays and updates the gas 
state. For example, the ray-drawing procedure can range in 
complexity from grids of rays restricted to radial paths origi- 
nating at the center of a spherical grid ( [Whalen and Norman, 
|2006| l, up to a variety of schemes for handling casting rays ei- 
ther with a fixed ( Abel et fl/.''1999'i 'Garcia-Segura and Franco 



middle ground ('Kessel-Deynet and Burkert 2000 Mac Low 

2006| l, while the most complex 



option is to restrict the hydrodynamic time step to resolve 



1996 ; Mac Low et al. , 2007 ; Melle ma ef a/. , 2006 ) or adaptive 
]Abel and Wandelt 2002; Krumh olz et al\ |2007c ; Rijkhorst 



et al. 2005j) r ay grid, or through a field of SPH particles ( ,Dale| 



TTal] |2007c| [Kessel-Deynet and Burkert[ [2000). Similarly 



there are a variety of possible time-stepping strategies for 
handling the interaction of radiation heating with (magneto- 
)hydrodynamics. The simplest are Stromgren volume meth- 
ods, in which one assumes that the gas reaches radiation 



and thermal equilibrium instantaneously (Dale et al. 2007c 



[Garcia-Segura and Franco[ [1996| l. Solving time-dependent 
equations for the thermal and chemical structure but not re- 
solving the relevant timescales hydrodynamically represents a 



gas heating and cooling times (Abel et al. 1999 Krumholz 
et a/.| p007a; Whalen and Norman 2006L As always in nu- 
merics, there is a tradeoff between speed and quality of so- 
lution. Fully resolving the ionization heating time produces 



measurably more accurate solutions ( Krumholz et al. 2007c I 



but is of course significantly more expensive than resolving it 
only marginally or assuming instantaneous equilibration. 



3. Numerical IVIethods for Protostellar Outflows 

Protostellar outflows are a natural result of the process of 
collapse and accretion that produces stars, and simulations of 
star formation that treat MHD and gravity with sufficient reso- 
lution do not need to include any additional physics to produce 
outflows (,Banerjee and Pudritz 2006 ; Machida et a/.[ [2008] ). 
However, sufficient resolution here means that the simulation 
must resolve the outflow launching region, which is typically 
no more than ^ 10 stellar radii. Achieving such high resolu- 
tion is prohibitively expensive for simulations that span more 
than a tiny fraction of the total formation time of a star, let 
alone an entire star cluster, so the most common procedure 
in such simulations is similar to that used for radiative feed- 
back. Replace the collapsing region with a sink of some sort, 
and model an outflow emerging from that sink via a subgrid 
model that sits on top of whatever sink model the computation 
uses. 

Such a model for outflows must specify the amount of mass 
and momentum contained, as well as the angular distribution 
of these quantities. Of these quantities the momentum is the 
best constrained by observations, since it remains unchanged 
even as the outflow gas entrains the material it encounters in 
the protostellar core. The mass flux and the angular distribu- 
tion of the outflow are less well constrained, since these are 
altered as the outflow ages and interacts with its environment. 
The correct value to use in a given simulation probably de- 
pends on the length scale that simulation resolves, since the 
mass and opening angle both increase as ouflowing gas moves 
away from the star and interacts with its environment. Most 
simulations assume the standard value of 10% for the ratio of 
infalling to outflowing mass. A good approximation for the 
opening angle is to assume the momentum of the outflowing 
gas is distributed with a profile p cx l/(r sin O^, where r is 
the distance from the star and 9 is the angle relative to the 
star's axis of rotation (Matzner and McKee 1999). The open- 
ing angle adopted in numerical simulations, however, varies 
enormously all the way from 0° ( Banerjee et al. 2001) to 90° 
( [Nakamura and Li[[2057l l. 

Once one has chosen a physical model for the outflow mass, 
momentum, and angular distribution in a given simulation, 
there remains the question of how to implement it numeri- 
cally. As noted above, all stars contribute significant amounts 
of momentum, so any realistic approach must include con- 
tributions from any star that forms in a simulation. In grid 
codes this problem is generally straightforward; one simply 
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adds mass and momentum with the desired angular distribu- 
tion to the computational cells in or immediately around the 
sink region for each star ( Nakamura and Li||2007| l. In particle- 
based codes the problem is more complex, since one must take 
care to avoid artificial clumping due to the discrete nature of 
the particles, and to ensure the momentum deposition in the 
region surrounding the sink is not altered by numerical inter- 
penetration of the SPH particles. A variety of strategies are 
available to solve these problems (Dale and Bonnell 2008| l, 
but they are computationally expensive, which presents a po- 
tential problem for simulations with large numbers of sources. 

Once the subgrid model is in place, outflows are much eas- 
ier to simulate than radiation feedback, because outflow evolu- 
tion is governed solely by hydrodynamics or MHD plus grav- 
ity, the physical mechanisms that are akeady included in any 
code used to simulate star formation. Beyond those involved 
in computing the subgrid model and injecting the outflow, the 
only additional computational cost that outflows impose on 
a code comes from the fact that outflow velocities can reach 
hundreds of km s~^, significantly greater than the ^ 10 km 
s^^ turbulent or infall speeds typically found in simulations 
that omit outflows. The higher speeds require smaller simula- 
tion time steps, and a corresponding increase in computational 
cost. 



V. SUMMARY AND OUTLOOK 

In this review we presented a overview of the current state 
of numerical star-formation studies. We have restricted our- 
selves to the early phases of stellar birth, from the formation 
of molecular clouds through to the build-up of stars and star 
clusters in their interior. We have left out the problem of accre- 
tion disks and protostellar evolution and point to other reviews 
in this context ( |Hartmann|| 19981 |Palla et a/.||2002[|Stahler and| 
[Pafla 2005). 

We hope we have illustrated that the question of stellar birth 
in our Galaxy and elsewhere in the universe is far from be- 
ing solved. Instead the field is rapidly evolving and has gone 
through a significant transformation in the last few years. In 
numerical star formation studies, we notice a general trend 
away from solely considering isolated processes and phenom- 
ena towards a more integrated multi-scale and multi-physics 
approach in todays computer simulations. In part this is 
triggered by the growing awareness that many physical pro- 
cesses contribute more or less equally to the formation of 
stars, such that it is not possible to single out individual ef- 
fects. Reliable and quantitative predictions can only be made 
on the basis of taking all relevant physical phenomena into 
account. Another reason for this development is the tremen- 
dous increase in computational capability provided by the ad- 
vent of (relatively) easy-to-handle massively-parallel super- 
computers, coupled with new and more efficient numerical al- 
gorithms for these machines. 

If we examine the past and current state of the art, then it 
is evident that most studies so far have focussed on a small 
number of physical processes only. Typically, one had a sin- 
gle question in mind, such as what happens if we include 



one particular physical process? How does it affect the sys- 
tem? How does it modify possible equilibrium states? And 
how does it influence the dynamical evolution if we apply 
perturbations? The processes and phenomena about which 
these questions have been asked include hydrodynamics, tur- 
bulence, gravitational dynamics, magnetic fields, nonequilib- 
rium chemistry, and the interaction of radiation with matter, 
but typically only one or two of them have been included 
in any given simulation. More sophisticated approaches in- 
clude larger numbers of processes, but no simulation so far 
has considered all of them. The challenge in the past was 
mainly to do justice to the inherent multidimensionality of the 
considered problems. For example, stellar birth in turbulent 
interstellar gas clouds with highly complex spatial and kine- 
matical structure is an intrinsically three-dimensional prob- 
lem with one- or two-dimensional approaches at best provid- 
ing order-of-magnitude estimates. Including multiple physi- 
cal processes gave way before the challenge of simulating in 
three dimensions. 

This era is coming to an end. Many of today's most chal- 
lenging problems are multi-physics, in the sense that they 
require the combination of many (if not all) of the above- 
mentioned processes, and multi-scale, in the sense that un- 
resolvable microscopic processes can feed back onto macro- 
scopic scales. This is true not only for star formation studies, 
but applies to virtually all fields of modern astrophysics. For 
example, the coagulation of dust species to larger particles or 
the interaction of dust with the radiation field from the central 
stars will eventually feedback into the dynamical behavior of 
the gas in protostellar accretion disks and hence has severe 
consequences for the formation and mass growths of plane- 
tary systems. Similarly, star formation and baryonic feedback 
are crucial ingredients of understanding galaxy formation and 
evolution in cosmological models. In a realistic description 
of cosmic phenomena, one is faced with the highly non-linear 
coupling between quite different kinds of interactions on a va- 
riety of scales. Star formation is no exception. 

This is not only a challenge, it is also a chance, because 
it may open up new pathways to successful collaborations 
across astrophysical disciplines. It also reaches out to sci- 
entists in neighboring fields, such as applied mathematics or 
computer science. For example, only few groups around the 
world are able to fully benefit from the massively parallel 
computing architectures that are currently being developed. 
Peak performances with ^100 teraflops will only be attain- 
able on thousands of CPUs, sustained petaflop computing may 
require as many as 10^ CPUs. This asks for a completely new 
approach to parallel algorithm design, a field where modern 
computer science is far ahead of the schemes currently used 
in astrophysics and star-formation studies. Regular method- 
ological exchange with applied mathematicians and possibly 
numerical fluid dynamicists thus holds the promise of both 
transferring new methods into astrophysics and raising the 
awareness of mathematicians about numerical challenges in 
astrophysics. 

Towards the end, we want to speculate about a few of the 
what we think are the most interesting and pressing open prob- 
lems in modern star formation theory and where the current 
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advancements in computational power and algorithmic so- 
phistication are likely to have a major impact. 

What drives interstellar turbulence? Observations show 
that turbulence in molecular clouds is ubiquitous, and that, 
with the exception of the dense cores discussed above, it 
seems to follow a universal relationship between velocity dis- 
persion and size. Even extragalactic molecular clouds appear 
to obey similar scalings. There are no variations in the turbu- 
lent properties between GMCs with little and much star for- 
mation, which might seem to argue for galaxy-scale driving, 
but there is also no systematic variation in GMC properties 
within a galaxy or between galaxies, which would seem to 
argue that internal processes must be important. So what is 
the relative importance of internal and external forcing mech- 
anisms in driving ISM turbulence? Does the answer depend 
on the length scales that one examines, or on the place where 
one looks? 

How does the multi-phase nature of the ISM influence stel- 
lar birth? Star formation appears to follow fairly universal 
scaling laws in galaxies that range from mildly Hl-dominated 
(such as the Milky Way) to galaxies that are strongly H2- 
dominated (such as local starbursts). Does the presence or 
absence of a significant atomic phase play an important role 
in regulating star formation, either directly (e.g. by limiting 
the amount of molecular gas "eligible" for star formation) or 
indirectly (e.g. by driving turbulent motions via thermal insta- 
bility)? How does the star formation process change, if at all, 
in galaxies such as dwarfs that contain very little molecular 
gas? 

How does stellar feedback influence star formation? Stars 
produce a wide variety of feedbacks: outflows, main sequence 
winds, ionizing and non-ionizing radiation, and supernovae. 
Which, if any of these, are responsible for controlling the 
rate and efficiency of star formation? Does the answer to 
this question change in different galactic environments, i.e. 
are there different processes acting in the denser molecular 
clouds found in circum-nuclear starbursts than in the tenuous 
outer regions of the galaxy? 

What determines the statistical properties of a stellar pop- 
ulation, and are these properties universal? On the observa- 
tional side, is the stellar IMF and binary distribution at present 
days different in different galactic environments, or is it truly 
universal? Especially in rich clusters our observational ba- 
sis still needs to be extended. The same holds for variations 
with metallicity as can be traced in the Local Group. Is the 
IMF in the Large Magellanic Cloud (with metal abundances 
of ~ 1/3 of the solar value) and the Small Magellanic cloud 
(with ~ 1/10 of that value) really similar to the Milky Way? 
On the theoretical side, what processes are responsible for 
the (non-)variation of the IMF? The critical mass for gravi- 
tational collapse can vary enormously between different en- 
vironments. Yet the IMF in globular clusters, for example, 
appears to be the same as in regions of distributed star for- 
mation as in Taurus. Hence, there must be additional physi- 
cal processes that influence the fragmentation behavior of the 
interstellar gas and determine the resulting stellar mass spec- 
trum. 
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